Connection check is rewritten to include the pi anchor bonds (may be incomplete yet)

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

diff --git a/MolLib/MD/MDCore.c b/MolLib/MD/MDCore.c
index 74ad09c..b79fcb7 100644 (file)
--- a/MolLib/MD/MDCore.c
+++ b/MolLib/MD/MDCore.c
@@ -276,10 +276,11 @@ s_register_missing_parameters(Int **missing, Int *nmissing, Int type, Int t1, In
 /*  Check the bonded atoms and append to results if not already present */
 /*  results[] is terminated by -1, hence must be at least (natom+1) size  */
 static int
-s_check_bonded(Atom *ap, Int *results)
+s_check_bonded(Molecule *mol, Int idx, Int *results)
 {
        Int i, n, *ip;
        const Int *cp;
+       Atom *ap = ATOM_AT_INDEX(mol->atoms, idx);
        cp = AtomConnectData(&ap->connect);
        for (i = 0; i < ap->connect.count; i++, cp++) {
                n = *cp;
@@ -292,6 +293,19 @@ s_check_bonded(Atom *ap, Int *results)
                        *ip = -1;
                }
        }
+       if (mol->piconnects != NULL) {
+               for (i = mol->piconnects[idx]; i < mol->piconnects[idx + 1]; i++) {
+                       n = mol->piconnects[i];
+                       for (ip = results; *ip >= 0; ip++) {
+                               if (n == *ip)
+                                       break;
+                       }
+                       if (*ip < 0) {
+                               *ip++ = n;
+                               *ip = -1;
+                       }
+               }
+       }
        for (ip = results; *ip >= 0; ip++)
                ;
        return ip - results;
@@ -302,7 +316,6 @@ s_make_exclusion_list(MDArena *arena)
 {
        Int *results;
        Int natoms = arena->mol->natoms;
-       Atom *atoms = arena->mol->atoms;
        MDExclusion *exinfo;
        int next_index, i, j;
 
@@ -332,15 +345,15 @@ s_make_exclusion_list(MDArena *arena)
                results[1] = -1;
                exinfo[i].index1 = 1;
                /*  1-2 exclusion (directly bonded)  */
-               exinfo[i].index2 = s_check_bonded(&atoms[i], results);
+               exinfo[i].index2 = s_check_bonded(arena->mol, i, results);
                n = exinfo[i].index2;
                /*  1-3 exclusion: atoms bonded to 1-2 exclusions  */
                for (j = exinfo[i].index1; j < exinfo[i].index2; j++)
-                       n = s_check_bonded(&atoms[results[j]], results);
+                       n = s_check_bonded(arena->mol, results[j], results);
                exinfo[i].index3 = n;
                /*  1-4 exclusion: atoms bonded to 1-3 exclusions  */
                for (j = exinfo[i].index2; j < exinfo[i].index3; j++)
-                       n = s_check_bonded(&atoms[results[j]], results);
+                       n = s_check_bonded(arena->mol, results[j], results);
                AssignArray(&arena->exlist, &arena->nexlist, sizeof(Int), next_index + n, NULL);
                memcpy(arena->exlist + next_index, results, n * sizeof(Int));
                exinfo[i].index0 += next_index;
@@ -1208,7 +1221,7 @@ s_find_pibond_parameters(MDArena *arena)
                                        ptype = -1;
                                        break;
                                } else
-                                       types[j] = mol->piatoms[j].type;
+                                       types[j] = mol->piatoms[ii].type;
                        } else if (ii < mol->natoms) {
                                idx[j] = ii;
                                types[j] = ATOM_AT_INDEX(mol->atoms, ii)->type;

diff --git a/MolLib/MD/MDForce.c b/MolLib/MD/MDForce.c
index 6a349fc..5e1c716 100644 (file)
--- a/MolLib/MD/MDForce.c
+++ b/MolLib/MD/MDForce.c
@@ -1379,6 +1379,7 @@ calc_force(MDArena *arena)
        s_calc_angle_force(arena);
        s_calc_dihedral_force(arena);
        s_calc_improper_force(arena);
+       s_calc_pibond_force(arena);
        s_calc_nonbonded_force(arena);
        s_calc_auxiliary_force(arena);
 

diff --git a/MolLib/MolAction.c b/MolLib/MolAction.c
index baf1550..5ee2931 100644 (file)
--- a/MolLib/MolAction.c
+++ b/MolLib/MolAction.c
@@ -1547,6 +1547,7 @@ s_MolActionInsertPiBonds(Molecule *mol, MolAction *action, MolAction **actp)
                InsertArray(&mol->pibonds, &mol->npibonds, sizeof(Int) * 4, j, 1, ip1 + i * 4);
        }
        *actp = MolActionNew(gMolActionRemovePiBonds, ig);
+       MoleculeInvalidatePiConnectionTable(mol);
        return 0;
 }
 
@@ -1565,6 +1566,7 @@ s_MolActionRemovePiBonds(Molecule *mol, MolAction *action, MolAction **actp)
        }
        *actp = MolActionNew(gMolActionInsertPiBonds, ig, n * 4, ip1);
        free(ip1);
+       MoleculeInvalidatePiConnectionTable(mol);
        return 0;
 }
 
@@ -1612,6 +1614,7 @@ s_MolActionInsertOnePiAtom(Molecule *mol, MolAction *action, MolAction **actp, i
                                mol->pibonds[i]++;
                }
        }
+       MoleculeInvalidatePiConnectionTable(mol);
        return 0;
 }
 
@@ -1646,6 +1649,7 @@ s_MolActionRemoveOnePiAtom(Molecule *mol, MolAction *action, MolAction **actp)
                if (mol->pibonds[i] > ATOMS_MAX_NUMBER + idx)
                        mol->pibonds[i]--;
        }       
+       MoleculeInvalidatePiConnectionTable(mol);
        return 0;
 }
 

diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c
index f059788..cc2565f 100755 (executable)
--- a/MolLib/Molecule.c
+++ b/MolLib/Molecule.c
@@ -376,6 +376,10 @@ MoleculeInitWithMolecule(Molecule *mp2, Molecule *mp)
                NewArray(&(mp2->pibonds), &(mp2->npibonds), sizeof(Int) * 4, mp->npibonds);
                memmove(mp2->pibonds, mp->pibonds, sizeof(Int) * 4 * mp->npibonds);
        }
+       if (mp->npiconnects > 0) {
+               NewArray(&(mp2->piconnects), &(mp2->npiconnects), sizeof(Int), mp->npiconnects);
+               memmove(mp2->piconnects, mp->piconnects, sizeof(Int) * mp->npiconnects);
+       }
        
 /*     mp2->useFlexibleCell = mp->useFlexibleCell; */
        if (mp->nframe_cells > 0) {
@@ -535,6 +539,11 @@ MoleculeClear(Molecule *mp)
                mp->pibonds = NULL;
                mp->npibonds = 0;
        }
+       if (mp->piconnects != NULL) {
+               free(mp->piconnects);
+               mp->piconnects = NULL;
+               mp->npiconnects = 0;
+       }
        if (mp->selection != NULL) {
                IntGroupRelease(mp->selection);
                mp->selection = NULL;
@@ -6467,6 +6476,7 @@ MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indice
                return -2;
 
        /*  Create atoms, with avoiding duplicates  */
+       MoleculeInvalidatePiConnectionTable(mp);
        n0 = n1 = mp->natoms;
        table = (int *)malloc(sizeof(int) * n0);
        if (table == NULL)
@@ -6842,6 +6852,7 @@ MoleculeCreateAnAtom(Molecule *mp, const Atom *ap, int pos)
                        if (mp->impropers[i] >= pos)
                                mp->impropers[i]++;
                }
+               MoleculeInvalidatePiConnectionTable(mp);
        }
        mp->nframes = -1;  /*  Should be recalculated later  */
        MoleculeIncrementModifyCount(mp);
@@ -6872,6 +6883,7 @@ MoleculeMerge(Molecule *dst, Molecule *src, IntGroup *where, int resSeqOffset)
                return 1;  /*  Bad parameter  */
 
        __MoleculeLock(dst);
+       MoleculeInvalidatePiConnectionTable(dst);
        nsrc = src->natoms;
        ndst = dst->natoms;
        if (resSeqOffset < 0)
@@ -7149,6 +7161,7 @@ sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqO
                return 1;  /*  Bad parameter  */
 
        __MoleculeLock(src);
+       MoleculeInvalidatePiConnectionTable(src);
        
        nsrc = src->natoms;
        nsrcnew = nsrc - ndst;
@@ -8773,6 +8786,8 @@ sMoleculeReorder(Molecule *mp)
        free(newAtoms);
        free(old2new);
        free(apArray);
+       
+       MoleculeInvalidatePiConnectionTable(mp);
 
 }
 
@@ -8876,6 +8891,8 @@ MoleculeRenumberAtoms(Molecule *mp, const Int *new2old, Int *old2new_out, Int is
        for (i = 0; i < mp->natoms; i++)
                memmove(ATOM_AT_INDEX(mp->atoms, old2new[i]), ATOM_AT_INDEX(saveAtoms, i), gSizeOfAtomRecord);
        retval = 0;
+       
+       MoleculeInvalidatePiConnectionTable(mp);
 
        MoleculeIncrementModifyCount(mp);
        mp->needsMDRebuild = 1;
@@ -9538,18 +9555,26 @@ static void
 sMoleculeFragmentSub(Molecule *mp, int idx, IntGroup *result, IntGroup *exatoms)
 {
        Atom *ap;
-       Int i, *cp;
+       Int i, *cp, idx2;
        if (exatoms != NULL && IntGroupLookup(exatoms, idx, NULL))
                return;
        IntGroupAdd(result, idx, 1);
        ap = ATOM_AT_INDEX(mp->atoms, idx);
        cp = AtomConnectData(&ap->connect);
        for (i = 0; i < ap->connect.count; i++) {
-               int idx2 = cp[i];
+               idx2 = cp[i];
                if (IntGroupLookup(result, idx2, NULL))
                        continue;
                sMoleculeFragmentSub(mp, idx2, result, exatoms);
        }
+       if (mp->piconnects != NULL) {
+               for (i = mp->piconnects[idx]; i < mp->piconnects[idx + 1]; i++) {
+                       idx2 = mp->piconnects[i];
+                       if (IntGroupLookup(result, idx2, NULL))
+                               continue;
+                       sMoleculeFragmentSub(mp, idx2, result, exatoms);
+               }
+       }
 }
 
 /*  The molecular fragment (= interconnected atoms) containing the atom n1 and
@@ -9561,6 +9586,7 @@ MoleculeFragmentExcludingAtomGroup(Molecule *mp, int n1, IntGroup *exatoms)
        if (mp == NULL || mp->natoms == 0 || n1 < 0 || n1 >= mp->natoms)
                return NULL;
        result = IntGroupNew();
+       MoleculeValidatePiConnectionTable(mp);
        sMoleculeFragmentSub(mp, n1, result, exatoms);
        return result;
 }
@@ -9578,6 +9604,7 @@ MoleculeFragmentExcludingAtoms(Molecule *mp, int n1, int argc, int *argv)
        for (i = 0; i < argc; i++)
                IntGroupAdd(exatoms, argv[i], 1);
        result = IntGroupNew();
+       MoleculeValidatePiConnectionTable(mp);
        sMoleculeFragmentSub(mp, n1, result, exatoms);
        IntGroupRelease(exatoms);
        return result;
@@ -9595,6 +9622,7 @@ MoleculeFragmentWithAtomGroups(Molecule *mp, IntGroup *inatoms, IntGroup *exatom
                return NULL;
        IntGroupIteratorInit(inatoms, &iter);
        result = IntGroupNew();
+       MoleculeValidatePiConnectionTable(mp);
        while ((i = IntGroupIteratorNext(&iter)) >= 0) {
                sMoleculeFragmentSub(mp, i, result, exatoms);
        }
@@ -9610,10 +9638,11 @@ MoleculeFragmentWithAtomGroups(Molecule *mp, IntGroup *inatoms, IntGroup *exatom
 int
 MoleculeIsFragmentDetachable(Molecule *mp, IntGroup *group, int *n1, int *n2)
 {
-       Int i, i1, i2, j, k, bond_count, nval1, nval2, *cp;
+       Int i, i1, i2, j, k, k2, bond_count, nval1, nval2, *cp;
        Atom *ap;
        if (mp == NULL || mp->natoms == 0 || group == NULL)
                return 0;  /*  Invalid arguments  */
+       MoleculeValidatePiConnectionTable(mp);
        bond_count = 0;
        for (i = 0; (i1 = IntGroupGetStartPoint(group, i)) >= 0; i++) {
                i2 = IntGroupGetEndPoint(group, i);
@@ -9631,6 +9660,18 @@ MoleculeIsFragmentDetachable(Molecule *mp, IntGroup *group, int *n1, int *n2)
                                                return 0;  /*  Too many bonds  */
                                }
                        }
+                       if (mp->piconnects != NULL) {
+                               for (k = mp->piconnects[j]; k < mp->piconnects[j + 1]; k++) {
+                                       k2 = mp->piconnects[k];
+                                       if (IntGroupLookup(group, k2, NULL) == 0) {
+                                               bond_count++;
+                                               nval1 = j;
+                                               nval2 = k2;
+                                               if (bond_count > 1)
+                                                       return 0;  /*  Too many bonds  */
+                                       }
+                               }
+                       }
                }
        }
        if (bond_count == 1) {
@@ -10069,6 +10110,95 @@ MoleculeCalculatePiAtomPosition(Molecule *mol, int idx, Vector *vp)
        return idx;
 }
 
+int
+MoleculeValidatePiConnectionTable(Molecule *mol)
+{
+       Int pass, i, j, k, m;
+       
+       if (mol == NULL || mol->pibonds == NULL)
+               return -1; /*  No need to process */
+       if (mol->piconnects != NULL)
+               return 0;  /*  Already valid  */
+
+       /*  Allocate index table  */
+       NewArray(&mol->piconnects, &mol->npiconnects, sizeof(Int), mol->natoms + 1);
+       memset(mol->piconnects, 0, sizeof(Int) * (mol->natoms + 1));
+       
+       /*  Pass 1: count connections for each atom  */
+       /*  Pass 2: store connection info  */
+       for (pass = 0; pass < 2; pass++) {
+               Int n[2], *ip[2], c[2];
+               for (i = 0; i < mol->npibonds; i++) {
+                       AtomConnect *ac;
+                       if (mol->pibonds[i * 4 + 2] >= 0)
+                               continue;  /*  Skip angle or dihedral entries  */
+                       for (j = 0; j < 2; j++) {
+                               n[j] = mol->pibonds[i * 4 + j];
+                               if (n[j] >= ATOMS_MAX_NUMBER) {
+                                       ac = &(mol->piatoms[n[j] - ATOMS_MAX_NUMBER].connect);
+                                       ip[j] = AtomConnectData(ac);
+                                       c[j] = ac->count;
+                               } else if (n[j] >= 0 && n[j] < mol->natoms) {
+                                       ip[j] = &n[j];
+                                       c[j] = 1;
+                               } else break;
+                       }
+                       if (j < 2)
+                               continue;  /*  Ignore the invalid entry  */
+                       for (j = 0; j < c[0]; j++) {
+                               for (k = 0; k < c[1]; k++) {
+                                       Int a1 = ip[0][j];
+                                       Int a2 = ip[1][k];
+                                       if (pass == 0) {
+                                               /*  Count  */
+                                               mol->piconnects[a1]++;
+                                               mol->piconnects[a2]++;
+                                       } else {
+                                               /*  Store the entry (at the first empty slot)  */
+                                               for (m = mol->piconnects[a1]; m < mol->piconnects[a1 + 1]; m++) {
+                                                       if (mol->piconnects[m] == -1) {
+                                                               mol->piconnects[m] = a2;
+                                                               break;
+                                                       }
+                                               }
+                                               for (m = mol->piconnects[a2]; m < mol->piconnects[a2 + 1]; m++) {
+                                                       if (mol->piconnects[m] == -1) {
+                                                               mol->piconnects[m] = a1;
+                                                               break;
+                                                       }
+                                               }
+                                       }
+                               }
+                       }
+               }
+               if (pass == 0) {
+                       /*  Expand the table, and store the position numbers  */
+                       m = mol->natoms + 1;
+                       for (i = 0; i <= mol->natoms; i++) {
+                               j = mol->piconnects[i];
+                               mol->piconnects[i] = m;
+                               m += j;
+                       }
+                       AssignArray(&mol->piconnects, &mol->npiconnects, sizeof(Int), m - 1, NULL);
+                       for (j = mol->natoms + 1; j < m; j++)
+                               mol->piconnects[j] = -1;
+               }
+       }
+       return (mol->npiconnects - mol->natoms - 1);  /*  Returns the number of entries  */
+}
+
+void
+MoleculeInvalidatePiConnectionTable(Molecule *mol)
+{
+       if (mol == NULL)
+               return;
+       if (mol->piconnects != NULL) {
+               free(mol->piconnects);
+               mol->piconnects = NULL;
+       }
+       mol->npiconnects = 0;
+}
+
 #pragma mark ====== MO calculation ======
 
 /*  Calculate an MO value for a single point.  */

diff --git a/MolLib/Molecule.h b/MolLib/Molecule.h
index 7f49250..579ba20 100755 (executable)
--- a/MolLib/Molecule.h
+++ b/MolLib/Molecule.h
@@ -277,7 +277,12 @@ typedef struct Molecule {
                                                    where n# is atom index if it is <ATOMS_MAX_NUMBER and
                                                    is piatom index + ATOMS_MAX_NUMBER otherwise.
                                                    The size of array is 4*npibonds.  */
-       
+       Int    npiconnects;  /*  Connection table for pi-metal bonds  */
+       Int    *piconnects;  /*  The first (natoms + 1) entries are for lookup table, and
+                                                    the connection data (indices of connected atoms) follow.
+                                                    The connected atoms (only by pi-bonds) for atom i are listed in
+                                                    piconnects[piconnects[i]]...piconnects[piconnects[i+1]]  */
+
        IntGroup *selection;
        Int    nframes;      /*  The number of frames (>= 1). This is a cached value, and should be
                                                         recalculated from the atoms if it is -1  */
@@ -514,7 +519,9 @@ int MoleculeSelectFrame(Molecule *mp, int frame, int copyback);
 int MoleculeFlushFrames(Molecule *mp);
 
 int MoleculeCalculatePiAtomPosition(Molecule *mol, int idx, Vector *vp);
-
+int MoleculeValidatePiConnectionTable(Molecule *mol);
+void MoleculeInvalidatePiConnectionTable(Molecule *mol);
+       
 int MoleculeCalcMO(Molecule *mp, Int mono, const Vector *op, const Vector *dxp, const Vector *dyp, const Vector *dzp, Int nx, Int ny, Int nz, int (*callback)(double progress, void *ref), void *ref);
 int MoleculeGetDefaultMOGrid(Molecule *mp, Int npoints, Vector *op, Vector *xp, Vector *yp, Vector *zp, Int *nx, Int *ny, Int *nz);
 const Cube *MoleculeGetCubeAtIndex(Molecule *mp, Int index);

diff --git a/MolLib/Ruby_bind/ruby_bind.c b/MolLib/Ruby_bind/ruby_bind.c
index f5fda03..8ff8576 100644 (file)
--- a/MolLib/Ruby_bind/ruby_bind.c
+++ b/MolLib/Ruby_bind/ruby_bind.c
@@ -4249,7 +4249,7 @@ s_Molecule_AtomOrPiAtomIndexFromValue(Molecule *mol, VALUE val)
                return n;
        p = StringValuePtr(val);
        for (i = 0; i < mol->npiatoms; i++) {
-               if (strcmp(p, mol->piatoms->aname) == 0)
+               if (strcmp(p, mol->piatoms[i].aname) == 0)
                        return ATOMS_MAX_NUMBER + i;  /*  PiAtom; should be looked up after normal atoms  */
        }
        if (n == -1)

diff --git a/Version b/Version
index afd04da..e0b9a3e 100644 (file)
--- a/Version
+++ b/Version
@@ -1,2 +1,2 @@
 version = "0.6.4"
-date = "20121023"
+date = "20121024"