Handling of periodic box during MD is being reworked.

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

diff --git a/MolLib/MD/MDCore.c b/MolLib/MD/MDCore.c
index 10f6c37..34ce958 100644 (file)
--- a/MolLib/MD/MDCore.c
+++ b/MolLib/MD/MDCore.c
@@ -1515,21 +1515,7 @@ md_prepare(MDArena *arena, int check_only)
 /*     if (t1 > mol->natoms && mol->box == NULL)
                return "symmetry operation is used but no periodic box is defined"; */
 
-       if (mol->cell != NULL) {
-               arena->periodic_a = (mol->cell->flags[0] != 0);
-               arena->periodic_b = (mol->cell->flags[1] != 0);
-               arena->periodic_c = (mol->cell->flags[2] != 0);
-               if (mol->nsyms > 0) {
-                       arena->cellsyms = (Transform *)realloc(arena->cellsyms, sizeof(Transform) * mol->nsyms);
-                       for (i = 0; i < mol->nsyms; i++) {
-                               Transform temp;
-                               TransformMul(temp, mol->syms[i], mol->cell->rtr);
-                               TransformMul(arena->cellsyms[i], mol->cell->tr, temp);
-                       }
-               }
-       } else {
-               arena->periodic_a = arena->periodic_b = arena->periodic_c = 0;
-       }
+       md_set_cell(arena);
 
        arena->natoms_uniq = t1;
        
@@ -2603,7 +2589,7 @@ md_minimize_step(MDArena *arena)
 }
 
 void
-md_cell_recalculate(MDArena *arena)
+md_update_cell(MDArena *arena)
 {
        Molecule *mol = arena->mol;
        XtalCell *cell = mol->cell;
@@ -2616,6 +2602,38 @@ md_cell_recalculate(MDArena *arena)
        }
 }
 
+void
+md_set_cell(MDArena *arena)
+{
+       Molecule *mol = arena->mol;
+       if (mol == NULL)
+               return;
+       if (mol->cell != NULL) {
+               arena->periodic_a = (mol->cell->flags[0] != 0);
+               arena->periodic_b = (mol->cell->flags[1] != 0);
+               arena->periodic_c = (mol->cell->flags[2] != 0);
+               if (mol->nsyms > 0) {
+                       Int i;
+                       if (mol->nsyms != arena->ncellsyms) {
+                               arena->ncellsyms = mol->nsyms;
+                               arena->cellsyms = (Transform *)realloc(arena->cellsyms, sizeof(Transform) * mol->nsyms);
+                       }
+                       for (i = 0; i < mol->nsyms; i++) {
+                               Transform temp;
+                               TransformMul(temp, mol->syms[i], mol->cell->rtr);
+                               TransformMul(arena->cellsyms[i], mol->cell->tr, temp);
+                       }
+               } else {
+                       arena->ncellsyms = 0;
+                       if (arena->cellsyms != NULL)
+                               free(arena->cellsyms);
+                       arena->cellsyms = NULL;
+               }
+       } else {
+               arena->periodic_a = arena->periodic_b = arena->periodic_c = 0;
+       }
+}
+
 /*  scale_atoms = 1: symmetry-unique atoms are transformed to keep the fractional coordinate constant */
 /*  scale_atoms = 2: same as scale_atoms = 1, except that the center of mass of each fragment of connected
        atoms are scaled and the atoms within one fragment are moved by the same amount */
@@ -2635,7 +2653,7 @@ md_scale_cell(MDArena *arena, const Transform tf, int scale_atoms)
        TransformVec(&cell->axes[0], tf, &cell->axes[0]);
        TransformVec(&cell->axes[1], tf, &cell->axes[1]);
        TransformVec(&cell->axes[2], tf, &cell->axes[2]);
-       md_cell_recalculate(arena);
+       md_update_cell(arena);
 
        /*  Deformation gradient (temp) = celltr * old_rcelltr  */
        TransformMul(grad, cell->tr, grad);
@@ -3002,7 +3020,7 @@ md_set_default(MDArena *arena)
        arena->cellc.y = 0;
        arena->cellc.z = 1; */
 /*     if (arena->mol != NULL && arena->mol->box != NULL)
-               md_cell_recalculate(arena); */
+               md_update_cell(arena); */
 }
 
 MDArena *

diff --git a/MolLib/MD/MDCore.h b/MolLib/MD/MDCore.h
index c55e9c5..19ebc35 100644 (file)
--- a/MolLib/MD/MDCore.h
+++ b/MolLib/MD/MDCore.h
@@ -335,6 +335,7 @@ typedef struct MDArena {
        /*  Symmetry operation  */
 //     Int    nsyms;
 //     Transform  *syms;
+       Int ncellsyms;
        Transform  *cellsyms;  /*  celltr * syms[i] * rcelltr  */
        
        /*  Unit cell vectors  */
@@ -452,7 +453,8 @@ int md_check_abnormal_improper(MDArena *arena, Molecule *mol, int idx);
 int md_set_alchemical_flags(MDArena *arena, int nflags, const char *flags);
 
 void md_amend_by_symmetry(MDArena *arena);
-void md_cell_recalculate(MDArena *arena);
+void md_set_cell(MDArena *arena);
+void md_update_cell(MDArena *arena);
 void md_scale_cell(MDArena *arena, const Transform tf, int scale_atoms);
 void md_wrap_coordinates(MDArena *arena);
 

diff --git a/MolLib/MD/MDPressure.c b/MolLib/MD/MDPressure.c
index 98f757c..cc4a438 100644 (file)
--- a/MolLib/MD/MDPressure.c
+++ b/MolLib/MD/MDPressure.c
@@ -231,7 +231,7 @@ pressure_control(MDArena *arena)
                        needs_cell_recalculate = 1;
                }
                if (needs_cell_recalculate)
-                       md_cell_recalculate(arena);
+                       md_update_cell(arena);
 
                /*  Recalc the energies  */
                md_amend_by_symmetry(arena);
@@ -250,7 +250,7 @@ pressure_control(MDArena *arena)
                        cell->axes[1] = cellb_save;
                        cell->axes[2] = cellc_save;
                        cell->origin = cello_save;
-                       md_cell_recalculate(arena);
+                       md_update_cell(arena);
                        if (pressure->control_algorithm != 1) {
                                md_restore(arena, 0);
                                arena->total_energy = total_energy_save;

diff --git a/MolLib/MolAction.c b/MolLib/MolAction.c
index c4e9cc8..e55314e 100644 (file)
--- a/MolLib/MolAction.c
+++ b/MolLib/MolAction.c
@@ -1581,15 +1581,21 @@ MolActionPerform(Molecule *mol, MolAction *action)
        } else if (strcmp(action->name, gMolActionSetCell) == 0) {
                if ((result = s_MolActionSetCell(mol, action, &act2)) != 0)
                        return result;
+               if (mol->arena != NULL)
+                       md_set_cell(mol->arena);
                needsRebuildMDArena = 1;
        } else if (strcmp(action->name, gMolActionSetBox) == 0) {
                if ((result = s_MolActionSetBox(mol, action, &act2)) != 0)
                        return result;
-               needsRebuildMDArena = 1;
+               if (mol->arena != NULL)
+                       md_set_cell(mol->arena);
+               needsSymmetryAmendment = 1;
        } else if (strcmp(action->name, gMolActionClearBox) == 0) {
                if ((result = s_MolActionSetBox(mol, NULL, &act2)) != 0)
                        return result;
-               needsRebuildMDArena = 1;
+               if (mol->arena != NULL)
+                       md_set_cell(mol->arena);
+               needsSymmetryAmendment = 1;
        } else if (strcmp(action->name, gMolActionEnableCellFlexibility) == 0) {
                if ((result = s_MolActionSetCellFlexibility(mol, action, &act2)) != 0)
                        return result;