# Conformer searching!

int i, j, k, l, n, nconformers;
int nsteps = option("Steps per Torsion","intspin",2,100,3,1);
int excludeh = option("Exclude X-X-X-H", "checkbox", 1);
int excludering = option("Exclude Rings", "checkbox", 1);
int minimise = option("Minimisation", "intcombo", "None,'Steepest Descent','Conjugate Gradient',MOPAC", 1)
double oldtors, delta = 360.0 / nsteps;

Model conf, m = aten.frame;
createPatterns();

if (m.zmatrix.ntorsions == 0)
{
	printf("No torsions in model. Nothing to do!\n");
	return;
}

printf("Model contains %i torsions.\n", m.zmatrix.ntorsions);

# Construct list of variable torsions
ZMatrixElement torsions[m.zmatrix.ntorsions];
int ntorsions = 0;
for (ZMatrixElement el = m.zmatrix.elements; el; ++el) if (el.torsionatom)
{
	# Exclude torsions for which either terminal atom is a hydrogen?
	if (excludeh && ((el.targetatom.z == 1) || (el.torsionatom.z == 1)))
	{
		printf("Torsion %i-%i-%i-%i excluded because it has terminal hydrogens...\n", el.targetatom.id, el.distanceatom.id, el.angleatom.id, el.torsionatom.id);
		continue;
	}
	# Exclude torsions present in cyclic systems
	#if 
	++ntorsions;
	torsions[ntorsions] = el;
}

printf("Total number of variable torsions is %i.\n", ntorsions);
nconformers = nsteps;
for (i=1; i<ntorsions; ++i) nconformers *= nsteps;
printf("Total number of conformers to be generated (with nsteps=%i) is %i.\n", nsteps, nconformers);

# Generate conformers.
# We have a list of variable torsions from a zmatrix of the original model, but changing zmatrix values does not
# rotate other attached atoms on 'j' and 'k' of the torsional zmatrixelement. We cannot change the original model
# with 'settorsion' since this will invalidate the zmatrix and do something nasty to our stored torsion pointers.
# So, in the loop we will copy the last generated configuration (or the original, if it is the first step) and
# paste it to a new model. We will then use settorsion() to modify this new model using the (equivalent) atom IDs
# from the stored zmtrixelement.
int counters[ntorsions] = 0;
m.selectAll();
m.copy();
for (int t = 1; t<=nconformers; ++t)
{
	# Create a new model and paste the contents
	Model conf = newModel(toa("Conformer%i",t));
	conf.paste();

	# Always rotate first torsion (this is the 'most quickly varying index')
	i = torsions[1].targetAtom.id;
	j = torsions[1].distanceAtom.id;
	k = torsions[1].angleAtom.id;
	l = torsions[1].torsionAtom.id;
	# Get current torsion value
	oldtors = torsion(i,j,k,l);
	setTorsion(i,j,k,l,oldtors+delta);
	++counters[1];

	# Increase other counters and rotate other torsions if necessary
	for (n=1; n<ntorsions; ++n)
	{
		if (counters[n] != nsteps) break;

		# Reset counter and rotate next torsion
		counters[n] = 0;
		i = torsions[n+1].targetAtom.id;
		j = torsions[n+1].distanceAtom.id;
		k = torsions[n+1].angleAtom.id;
		l = torsions[n+1].torsionAtom.id;

		# Get current torsion value
		oldtors = torsion(i,j,k,l);
		setTorsion(i,j,k,l,oldtors+delta);

		# Minimise?
		if (minimise == 2) sdMinimise();
		else if (minimise == 3) cgMinimise();
		else if (minimise == 4) mopacMinimise();

		++counters[n+1];
	}

	#printf("COUNTERS = %i %i %i \n", counters[1], counters[2], counters[3]);

	# Save files...

}
