Cell minimization was not working correctly when no symmetry operations are defined...

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

diff --git a/MolLib/MD/MDCore.c b/MolLib/MD/MDCore.c
index a064910..6512f88 100644 (file)
--- a/MolLib/MD/MDCore.c
+++ b/MolLib/MD/MDCore.c
@@ -2284,9 +2284,9 @@ md_wrap_coordinates(MDArena *arena)
                        /*  Calculate the offset from the position of the center of mass */
                        Vector r = arena->fragment_info[n].pos;
                        TransformVec(&r, arena->mol->cell->rtr, &r);
-                       if (ap->symop.alive && ap->symop.sym > 0 && ap->symop.sym <= arena->mol->nsyms) {
+                       if (ap->symop.alive && ap->symop.sym > 0 && ap->symop.sym < arena->mol->nsyms) {
                                /*  The translational components of symop are not included  */
-                               TransformVec(&r, arena->mol->syms[ap->symop.sym - 1], &r);
+                               TransformVec(&r, SYMMETRY_AT_INDEX(arena->mol->syms, ap->symop.sym - 1), &r);
                        }
                        ap->wrap_dx = (arena->periodic_a ? -floor(r.x) : 0);
                        ap->wrap_dy = (arena->periodic_b ? -floor(r.y) : 0);
@@ -2873,7 +2873,7 @@ s_md_modify_cell_parameters(MDArena *arena, Double lambda)
                cell->axes[1].y = arena->old_cell_pars[4] + arena->cell_vels[4] * lambda;
                cell->axes[1].z = arena->old_cell_pars[5] + arena->cell_vels[5] * lambda;
        }
-       if (arena->periodic_b) {
+       if (arena->periodic_c) {
                cell->axes[2].x = arena->old_cell_pars[6] + arena->cell_vels[6] * lambda;
                cell->axes[2].y = arena->old_cell_pars[7] + arena->cell_vels[7] * lambda;
                cell->axes[2].z = arena->old_cell_pars[8] + arena->cell_vels[8] * lambda;
@@ -2945,8 +2945,8 @@ md_minimize_cell_step(MDArena *arena)
                return 1;  /*  Gradient is sufficiently small  */
        
        /*  Proceed along cell_vels[] until the energy increases  */
-       low_limit = arena->coordinate_convergence / arena->max_gradient;
-       high_limit = 0.1 / arena->max_gradient;
+       low_limit = arena->coordinate_convergence / arena->cell_max_gradient;
+       high_limit = 0.1 / arena->cell_max_gradient;
        low = 0.0;
        low_energy = arena->total_energy;
        lambda = high_limit;

diff --git a/MolLib/MD/MDForce.c b/MolLib/MD/MDForce.c
index 3fd87a6..f711b2a 100644 (file)
--- a/MolLib/MD/MDForce.c
+++ b/MolLib/MD/MDForce.c
@@ -936,7 +936,7 @@ s_calc_nonbonded_force_sub(MDArena *arena, Double *energies, Double *eenergies,
                                rij2 = ATOM_AT_INDEX(arena->mol->atoms, ap2->symbase)->r;
                        else rij2 = ap2->r;
                        TransformVec(&rr, cell->rtr, &rij2);
-                       memmove(symtr, arena->mol->syms[ap2->symop.sym], sizeof(Transform));
+                       memmove(symtr, SYMMETRY_AT_INDEX(arena->mol->syms, ap2->symop.sym), sizeof(Transform));
                        symtr[9] += ap2->symop.dx + vl->symop.dx;
                        symtr[10] += ap2->symop.dy + vl->symop.dy;
                        symtr[11] += ap2->symop.dz + vl->symop.dz;

diff --git a/MolLib/MainView.c b/MolLib/MainView.c
index fdd518f..9625081 100755 (executable)
--- a/MolLib/MainView.c
+++ b/MolLib/MainView.c
@@ -1207,7 +1207,7 @@ drawAtom(MainView *mview, int i1, int selected, const Vector *dragOffset, const
                ap = ATOM_AT_INDEX(mview->mol->atoms, ep->index);
                an1 = ap->atomicNumber;
                r1 = ap->r;
-               trp = &(mview->mol->syms[ep->symop]);
+               trp = &(SYMMETRY_AT_INDEX(mview->mol->syms, ep->symop));
                if (/* !mview->mol->is_xtal_coord && */ mview->mol->cell != NULL) {
                        TransformVec(&r1, mview->mol->cell->rtr, &r1);
                        TransformVec(&r1, *trp, &r1);
@@ -1346,7 +1346,7 @@ drawBond(MainView *mview, int i1, int i2, int selected, int selected2, int draft
                        ap = ATOM_AT_INDEX(mview->mol->atoms, ep->index);
                        an[i] = ap->atomicNumber;
                        r[i] = ap->r;
-                       TransformVec(&r[i], mview->mol->syms[ep->symop], &r[i]);
+                       TransformVec(&r[i], SYMMETRY_AT_INDEX(mview->mol->syms, ep->symop), &r[i]);
                        VecInc(r[i], ep->dr);
                } else {
                        ap = ATOM_AT_INDEX(mview->mol->atoms, in);

diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c
index 87b6126..72aa54c 100755 (executable)
--- a/MolLib/Molecule.c
+++ b/MolLib/Molecule.c
@@ -358,8 +358,7 @@ MoleculeInitWithMolecule(Molecule *mp2, const Molecule *mp)
                memmove(mp2->cell, mp->cell, sizeof(XtalCell));
        }
        if (mp->nsyms > 0) {
-               mp2->nsyms = mp->nsyms;
-               mp2->syms = (Transform *)calloc(sizeof(Transform), mp2->nsyms);
+               NewArray(&(mp2->syms), &(mp2->nsyms), sizeof(Transform), mp->nsyms);
                memmove(mp2->syms, mp->syms, sizeof(Transform) * mp2->nsyms);
        }
        mp2->useFlexibleCell = mp->useFlexibleCell;
@@ -6222,9 +6221,9 @@ MoleculeGetTransformForSymop(Molecule *mp, Symop symop, Transform *tf, int is_ca
        Transform t;
        if (mp == NULL || mp->cell == NULL)
                return -1;
-       if (symop.sym >= mp->nsyms)
+       if (symop.sym >= mp->nsyms && symop.sym != 0)
                return -2;
-       memmove(*tf, mp->syms[symop.sym], sizeof(Transform));
+       memmove(*tf, SYMMETRY_AT_INDEX(mp->syms, symop.sym), sizeof(Transform));
        (*tf)[9] += symop.dx;
        (*tf)[10] += symop.dy;
        (*tf)[11] += symop.dz;
@@ -6283,17 +6282,17 @@ MoleculeTransformBySymop(Molecule *mp, const Vector *vpin, Vector *vpout, Symop
 {
        if (mp == NULL)
                return 1;
-       if (symop.sym >= mp->nsyms)
+       if (symop.sym >= mp->nsyms && symop.sym != 0)
                return 2;
        if (mp->cell != NULL /* && !mp->is_xtal_coord */) {
                TransformVec(vpout, mp->cell->rtr, vpin);
-               TransformVec(vpout, mp->syms[symop.sym], vpout);
+               TransformVec(vpout, SYMMETRY_AT_INDEX(mp->syms, symop.sym), vpout);
                vpout->x += symop.dx;
                vpout->y += symop.dy;
                vpout->z += symop.dz;
                TransformVec(vpout, mp->cell->tr, vpout);
        } else {
-               TransformVec(vpout, mp->syms[symop.sym], vpin);
+               TransformVec(vpout, SYMMETRY_AT_INDEX(mp->syms, symop.sym), vpin);
                vpout->x += symop.dx;
                vpout->y += symop.dy;
                vpout->z += symop.dz;
@@ -9157,7 +9156,7 @@ MoleculeSetAnisoBySymop(Molecule *mp, int idx)
                memmove(ap->aniso, ap2->aniso, sizeof(Aniso));
                return;
        }
-       memmove(t1, mp->syms[ap->symop.sym], sizeof(Transform));
+       memmove(t1, SYMMETRY_AT_INDEX(mp->syms, ap->symop.sym), sizeof(Transform));
        t1[9] = t1[10] = t1[11] = 0.0;
        memset(t2, 0, sizeof(Transform));
        t2[0] = ap2->aniso->bij[0];

diff --git a/MolLib/Molecule.h b/MolLib/Molecule.h
index bf8482b..35aa31c 100755 (executable)
--- a/MolLib/Molecule.h
+++ b/MolLib/Molecule.h
@@ -110,6 +110,7 @@ extern Int gSizeOfAtomRecord;
 #define ATOM_PREV(p)         ((Atom *)((char *)(p) - gSizeOfAtomRecord))
 #define SYMOP_ALIVE(s) ((s.dx || s.dy || s.dz || s.sym) != 0)
 #define SYMOP_EQUAL(s1, s2) (s1.dx == s2.dx && s1.dy == s2.dy && s1.dz == s2.dz && s1.sym == s2.sym)
+#define SYMMETRY_AT_INDEX(p, i) (*((i) == 0 ? &gIdentityTransform : &p[i]))
 
 /*  atom.connects is a union entry, including direct data for nconnects <= ATOM_CONNECTS_LIMIT
     and malloc()'ed entry for nconnects > ATOM_CONNECTS_LIMIT. The following functions

diff --git a/MolLib/Types.c b/MolLib/Types.c
index 45ec8d1..e7b041f 100755 (executable)
--- a/MolLib/Types.c
+++ b/MolLib/Types.c
@@ -341,6 +341,8 @@ TransformForRotation(Transform dst, const Vector *axis, Double angle, const Vect
        return 0;
 }
 
+Transform gIdentityTransform = {1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0};
+
 #pragma mark ==== LAMatrix ====
 
 typedef struct LAMatrixTempRecord {
@@ -610,7 +612,9 @@ NewArray(void *base, Int *count, int item_size, int nitems)
        void **bp = (void *)base;
        *bp = NULL;
        *count = 0;
-       return AssignArray(base, count, item_size, nitems - 1, NULL);
+       if (nitems > 0)
+               return AssignArray(base, count, item_size, nitems - 1, NULL);
+       else return NULL;
 }
 
 /*  Insert items to an array.  */

diff --git a/MolLib/Types.h b/MolLib/Types.h
index dd78283..3c36863 100755 (executable)
--- a/MolLib/Types.h
+++ b/MolLib/Types.h
@@ -123,6 +123,8 @@ void TransformForInversion(Transform dst, const Vector *center);
 int  TransformForReflection(Transform dst, const Vector *axis, const Vector *center);
 int  TransformForRotation(Transform dst, const Vector *axis, Double angle, const Vector *center);
 
+extern Transform gIdentityTransform;
+
 /*  Wrapper struct for CLAPACK routines  */
 typedef struct LAMatrix {
        __CLPK_integer row, column;