From: toshinagata1964 <toshinagata1964@a2be9bc6-48de-4e38-9406-05402d4bc13c>
Date: Fri, 15 Nov 2013 02:45:35 +0000 (+0000)
Subject: Minimizing structure including pi-anchor atoms was not working correctly. Fixed.
X-Git-Tag: v1.0.2~221
X-Git-Url: http://git.osdn.net/view?a=commitdiff_plain;h=e25785c3eb89e10946010fa89d32a61b8d3300c5;p=molby%2FMolby.git

Minimizing structure including pi-anchor atoms was not working correctly. Fixed.

git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@417 a2be9bc6-48de-4e38-9406-05402d4bc13c
---

diff --git a/MolLib/MD/MDCore.c b/MolLib/MD/MDCore.c
index 88bff3a..a6e5a43 100644
--- a/MolLib/MD/MDCore.c
+++ b/MolLib/MD/MDCore.c
@@ -2198,23 +2198,8 @@ md_update_positions(MDArena *arena)
 	}
 
 	/*  Update the anchor positions  */
-	if (count_anchors > 0) {
-		Int *cp, n, j;
-		Atom *ap2;
-		ap = arena->mol->atoms;
-		for (i = 0; i < natoms; i++, ap++) {
-			if (ap->anchor == NULL)
-				continue;
-			cp = AtomConnectData(&ap->anchor->connect);
-			n = ap->anchor->connect.count;
-			VecZero(ap->r);
-			for (j = 0; j < n; j++) {
-				Double w = ap->anchor->coeffs[j];
-				ap2 = arena->mol->atoms + cp[j];
-				VecScaleInc(ap->r, ap2->r, w);
-			}
-		}
-	}
+	if (count_anchors > 0)
+		MoleculeUpdatePiAnchorPositions(arena->mol);
 	
 	/*  Update the abnormal bond parameters  */
 	if (arena->anbond_thres > 0.0) {
@@ -2505,7 +2490,7 @@ md_minimize_atoms_step(MDArena *arena)
 	Double bk, w1, w2, w3, dump;
 	Double low, mid, high, low_energy, mid_energy, high_energy, lambda;
 	Double low_limit, high_limit;
-	Int i, j, retval, natoms_movable;
+	Int i, j, retval, natoms_movable, nanchors;
 	Atom *atoms = arena->mol->atoms;
 	Atom *ap;
 	Int natoms = arena->mol->natoms;
@@ -2517,9 +2502,14 @@ md_minimize_atoms_step(MDArena *arena)
 	w1 = w2 = 0.0;
 	retval = 0;
 	natoms_movable = 0;
+	nanchors = 0;
 	for (i = 0, ap = atoms, vp = arena->old_forces; i < natoms; i++, ap++, vp++) {
 		if (ap->fix_force < 0)
 			continue;
+		if (ap->anchor != NULL) {
+			nanchors++;
+			continue;
+		}
 		w1 += VecLength2(ap->f);
 		w2 += VecDot(ap->f, *vp);
 		natoms_movable++;
@@ -2549,7 +2539,7 @@ md_minimize_atoms_step(MDArena *arena)
 	}
 	dump = sqrt(dump);
 	for (i = 0, ap = atoms, vp = arena->old_forces; i < natoms; i++, ap++, vp++) {
-		if (ap->fix_force < 0)
+		if (ap->fix_force < 0 || ap->anchor != NULL)
 			continue;
 		*vp = ap->f;
 		*(vp + natoms) = ap->r;
@@ -2563,6 +2553,8 @@ md_minimize_atoms_step(MDArena *arena)
 		if (w1 > w3)
 			w3 = w1;
 	}
+	if (nanchors > 0)
+		MoleculeUpdatePiAnchorPositions(arena->mol);
 	w3 = sqrt(w3);
 	arena->max_gradient = w3;
 	arena->v_len2 = w2;
@@ -2580,7 +2572,7 @@ md_minimize_atoms_step(MDArena *arena)
 	high = lambda;
 	while (1) {
 		for (j = 0, ap = atoms, vp = arena->old_pos, vdr = arena->verlets_dr; j < natoms; j++, ap++, vp++, vdr++) {
-			if (ap->fix_force < 0)
+			if (ap->fix_force < 0 || ap->anchor != NULL)
 				continue;
 			r = ap->r;
 			ap->r.x = vp->x + ap->v.x * lambda;
@@ -2589,6 +2581,8 @@ md_minimize_atoms_step(MDArena *arena)
 			VecDec(r, ap->r);
 			VecInc(*vdr, r);
 		}
+		if (nanchors > 0)
+			MoleculeUpdatePiAnchorPositions(arena->mol);
 		calc_force(arena);
 		mid = lambda;
 		mid_energy = arena->total_energy;
@@ -2608,13 +2602,15 @@ md_minimize_atoms_step(MDArena *arena)
 			/*  Cannot find point with lower energy than the starting point  */
 			/*  Restore the original position  */
 			for (j = 0, ap = atoms, vp = arena->old_pos, vdr = arena->verlets_dr; j < natoms; j++, ap++, vp++, vdr++) {
-				if (ap->fix_force < 0)
+				if (ap->fix_force < 0 || ap->anchor != NULL)
 					continue;
 				r = ap->r;
 				ap->r = *vp;
 				VecDec(r, ap->r);
 				VecInc(*vdr, r);
 			}
+			if (nanchors > 0)
+				MoleculeUpdatePiAnchorPositions(arena->mol);
 			calc_force(arena);
 			lambda = 0.0;
 			if (bk == 0.0)
@@ -2629,7 +2625,7 @@ md_minimize_atoms_step(MDArena *arena)
 		if (high - mid > mid - low) {
 			lambda = high - (high - mid) * phi;
 			for (j = 0, ap = atoms, vp = arena->old_pos, vdr = arena->verlets_dr; j < natoms; j++, ap++, vp++, vdr++) {
-				if (ap->fix_force < 0)
+				if (ap->fix_force < 0 || ap->anchor != NULL)
 					continue;
 				r = ap->r;
 				ap->r.x = vp->x + ap->v.x * lambda;
@@ -2638,6 +2634,8 @@ md_minimize_atoms_step(MDArena *arena)
 				VecDec(r, ap->r);
 				VecInc(*vdr, r);
 			}	
+			if (nanchors > 0)
+				MoleculeUpdatePiAnchorPositions(arena->mol);
 			calc_force(arena);
 			if (arena->total_energy < mid_energy) {
 				low = mid;
@@ -2651,7 +2649,7 @@ md_minimize_atoms_step(MDArena *arena)
 		} else {
 			lambda = mid - (mid - low) * phi;
 			for (j = 0, ap = atoms, vp = arena->old_pos, vdr = arena->verlets_dr; j < natoms; j++, ap++, vp++, vdr++) {
-				if (ap->fix_force < 0)
+				if (ap->fix_force < 0 || ap->anchor != NULL)
 					continue;
 				r = ap->r;
 				ap->r.x = vp->x + ap->v.x * lambda;
@@ -2660,6 +2658,8 @@ md_minimize_atoms_step(MDArena *arena)
 				VecDec(r, ap->r);
 				VecInc(*vdr, r);
 			}	
+			if (nanchors > 0)
+				MoleculeUpdatePiAnchorPositions(arena->mol);
 			calc_force(arena);
 			if (arena->total_energy < mid_energy) {
 				high = mid;
diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c
index 66c0bcf..53b8384 100755
--- a/MolLib/Molecule.c
+++ b/MolLib/Molecule.c
@@ -10473,23 +10473,43 @@ MoleculeFlushFrames(Molecule *mp)
 
 #pragma mark ====== Pi Atoms ======
 
+static inline void
+sMoleculeCalculatePiAnchorPosition(Atom *ap, Atom *atoms)
+{
+	Int *cp, j, n;
+	Atom *ap2;
+	cp = AtomConnectData(&ap->anchor->connect);
+	n = ap->anchor->connect.count;
+	VecZero(ap->r);
+	for (j = 0; j < n; j++) {
+		Double w = ap->anchor->coeffs[j];
+		ap2 = ATOM_AT_INDEX(atoms, cp[j]);
+		VecScaleInc(ap->r, ap2->r, w);
+	}	
+}
+
+void
+MoleculeUpdatePiAnchorPositions(Molecule *mol)
+{
+	Int i;
+	Atom *ap;
+	for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
+		if (ap->anchor == NULL)
+			continue;
+		sMoleculeCalculatePiAnchorPosition(ap, mol->atoms);
+	}
+}
+
 void
 MoleculeCalculatePiAnchorPosition(Molecule *mol, int idx)
 {
-	Atom *ap, *ap2;
-	Int i, n, *ip;
+	Atom *ap;
 	if (mol == NULL || idx < 0 || idx >= mol->natoms)
 		return;
 	ap = ATOM_AT_INDEX(mol->atoms, idx);
 	if (ap->anchor == NULL)
 		return;
-	ip = AtomConnectData(&ap->anchor->connect);
-	n = ap->anchor->connect.count;
-	VecZero(ap->r);
-	for (i = 0; i < ap->anchor->connect.count; i++) {
-		ap2 = ATOM_AT_INDEX(mol->atoms, ip[i]);
-		VecScaleInc(ap->r, ap2->r, ap->anchor->coeffs[i]);
-	}
+	sMoleculeCalculatePiAnchorPosition(ap, mol->atoms);
 }
 
 int
diff --git a/MolLib/Molecule.h b/MolLib/Molecule.h
index c895c91..a4468f2 100755
--- a/MolLib/Molecule.h
+++ b/MolLib/Molecule.h
@@ -522,6 +522,7 @@ int MoleculeRemoveFrames(Molecule *mp, IntGroup *group, Vector *outFrame, Vector
 int MoleculeSelectFrame(Molecule *mp, int frame, int copyback);
 int MoleculeFlushFrames(Molecule *mp);
 
+void MoleculeUpdatePiAnchorPositions(Molecule *mol);
 void MoleculeCalculatePiAnchorPosition(Molecule *mol, int idx);
 int MoleculeSetPiAnchorList(Molecule *mol, Int idx, Int nentries, Int *entries, Double *weights, Int *nUndoActions, struct MolAction ***undoActions);