Another implementation of pi-anchor is being tried (on the way; doesn't even compile...

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

diff --git a/MolLib/MD/MDCore.c b/MolLib/MD/MDCore.c
index b79fcb7..d5419ea 100644 (file)
--- a/MolLib/MD/MDCore.c
+++ b/MolLib/MD/MDCore.c
@@ -293,6 +293,7 @@ s_check_bonded(Molecule *mol, Int idx, Int *results)
                        *ip = -1;
                }
        }
+#if PIATOM
        if (mol->piconnects != NULL) {
                for (i = mol->piconnects[idx]; i < mol->piconnects[idx + 1]; i++) {
                        n = mol->piconnects[i];
@@ -306,6 +307,7 @@ s_check_bonded(Molecule *mol, Int idx, Int *results)
                        }
                }
        }
+#endif
        for (ip = results; *ip >= 0; ip++)
                ;
        return ip - results;
@@ -560,6 +562,7 @@ s_search_improper(MDArena *arena, int n1, int n2, int n3, int n4)
 
 #define SWAP_INT(_i, _j) do { Int _k; _k = _i; _i = _j; _j = _k; } while (0)
 
+#if PIATOM
 /*  Check whether i1-i2 is included in the pi-bond list  */
 static int
 s_check_pi_bond(MDArena *arena, Int i1, Int i2)
@@ -593,6 +596,7 @@ s_check_pi_bond(MDArena *arena, Int i1, Int i2)
        }
        return 0;
 }
+#endif
 
 /*  Find vdw parameters and build the in-use list */
 static int
@@ -832,11 +836,13 @@ s_find_bond_parameters(MDArena *arena)
                        Atom *ap1, *ap2;
                        i1 = mol->bonds[i * 2];
                        i2 = mol->bonds[i * 2 + 1];
+#if PIATOM
                        if (s_check_pi_bond(arena, i1, i2)) {
                                /*  Skip this bond  */
                                arena->bond_par_i[i] = -1;
                                continue;
                        }
+#endif
                        ap1 = ATOM_AT_INDEX(mol->atoms, i1);
                        ap2 = ATOM_AT_INDEX(mol->atoms, i2);
                        type1 = ap1->type;
@@ -924,11 +930,13 @@ s_find_angle_parameters(MDArena *arena)
                        i1 = mol->angles[i * 3];
                        i2 = mol->angles[i * 3 + 1];
                        i3 = mol->angles[i * 3 + 2];
+#if PIATOM
                        if (s_check_pi_bond(arena, i1, i2) || s_check_pi_bond(arena, i2, i3)) {
                                /*  Skip this angle  */
                                arena->angle_par_i[i] = -1;
                                continue;
                        }
+#endif
                        type1 = ATOM_AT_INDEX(mol->atoms, i1)->type;
                        type2 = ATOM_AT_INDEX(mol->atoms, i2)->type;
                        type3 = ATOM_AT_INDEX(mol->atoms, i3)->type;
@@ -1017,11 +1025,13 @@ s_find_dihedral_parameters(MDArena *arena)
                        i2 = mol->dihedrals[i * 4 + 1];
                        i3 = mol->dihedrals[i * 4 + 2];
                        i4 = mol->dihedrals[i * 4 + 3];
+#if PIATOM
                        if (s_check_pi_bond(arena, i1, i2) || s_check_pi_bond(arena, i2, i3) || s_check_pi_bond(arena, i3, i4)) {
                                /*  Skip this dihedral  */
                                arena->dihedral_par_i[i] = -1;
                                continue;
                        }
+#endif
                        type1 = ATOM_AT_INDEX(mol->atoms, i1)->type;
                        type2 = ATOM_AT_INDEX(mol->atoms, i2)->type;
                        type3 = ATOM_AT_INDEX(mol->atoms, i3)->type;
@@ -1188,6 +1198,7 @@ s_find_improper_parameters(MDArena *arena)
        return nmissing;
 }
 
+#if PIATOM
 /*  Find pi-bond parameters  */
 static int
 s_find_pibond_parameters(MDArena *arena)
@@ -1293,6 +1304,7 @@ s_find_pibond_parameters(MDArena *arena)
        arena->nmissing += nmissing;
        return nmissing;
 }
+#endif
 
 /*  Find one fragment, starting from start_index  */
 static void
@@ -1707,10 +1719,8 @@ md_prepare(MDArena *arena, int check_only)
                   "Number of bonds = %d\n"
                   "Number of angles = %d\n"
                   "Number of dihedrals = %d\n"
-                  "Number of impropers = %d\n"
-                  "Number of pi anchors = %d\n"
-                  "Number of pi anchor constructs = %d\n",
-                  mol->nbonds, mol->nangles, mol->ndihedrals, mol->nimpropers, mol->npiatoms, mol->npibonds);
+                  "Number of impropers = %d\n",
+                  mol->nbonds, mol->nangles, mol->ndihedrals, mol->nimpropers);
        
        t1 = t2 = t3 = 0;
        for (i = 0; i < mol->natoms; i++) {
@@ -1744,19 +1754,24 @@ md_prepare(MDArena *arena, int check_only)
        if (arena->par != NULL)
                ParameterRelease(arena->par);
        arena->par = ParameterNew();
+
+#if PIATOM
        if (arena->pi_pars != NULL)
                free(arena->pi_pars);
        if (arena->mol->npibonds > 0)
                arena->pi_pars = (UnionPar *)malloc(sizeof(UnionPar) * arena->mol->npibonds);
        else arena->pi_pars = NULL;
-
+#endif
+       
        arena->nmissing = arena->nsuspicious = 0;
        s_find_vdw_parameters(arena);
        s_find_bond_parameters(arena);
        s_find_angle_parameters(arena);
        s_find_dihedral_parameters(arena);
        s_find_improper_parameters(arena);
+#if PIATOM
        s_find_pibond_parameters(arena);
+#endif
        
        if (arena->nmissing > 0) {
        /*      for (i = 0; i < nmissing; i++) {
@@ -3725,12 +3740,14 @@ md_arena_set_molecule(MDArena *arena, Molecule *xmol)
                memmove(mol->impropers, xmol->impropers, sizeof(Int) * 4 * xmol->nimpropers);
                NewArray(&mol->syms, &mol->nsyms, sizeof(Transform), xmol->nsyms);
                memmove(mol->syms, xmol->syms, sizeof(Transform) * xmol->nsyms);
+#if PIATOM
                NewArray(&mol->piatoms, &mol->npiatoms, sizeof(PiAtom), xmol->npiatoms);
                for (i = 0; i < xmol->npiatoms; i++) {
                        PiAtomDuplicate(mol->piatoms + i, xmol->piatoms + i);
                }
                NewArray(&mol->pibonds, &mol->npibonds, sizeof(Int) * 4, xmol->npibonds);
                memmove(mol->pibonds, xmol->pibonds, sizeof(Int) * 4 * xmol->npibonds);
+#endif
                if (xmol->cell != NULL) {
                        mol->cell = (XtalCell *)malloc(sizeof(XtalCell));
                        memmove(mol->cell, xmol->cell, sizeof(XtalCell));

diff --git a/MolLib/MD/MDForce.c b/MolLib/MD/MDForce.c
index 31779fe..d7f6754 100644 (file)
--- a/MolLib/MD/MDForce.c
+++ b/MolLib/MD/MDForce.c
@@ -468,6 +468,7 @@ s_calc_improper_force(MDArena *arena)
        s_calc_dihedral_force_sub(arena, mol->atoms, mol->nimpropers, mol->impropers, arena->improper_par_i, par->improperPars, NULL/*arena->sym_improper_uniq*/, energies, forces);
 }
 
+#if PIATOM
 static void
 s_calc_pibond_force(MDArena *arena)
 {
@@ -552,6 +553,7 @@ s_calc_pibond_force(MDArena *arena)
                }
        }
 }
+#endif  /*  PIATOM  */
 
 /*  ==================================================================== */
 /*  md_check_verlet_list: only for debugging the verlet list generation  */
@@ -1382,7 +1384,9 @@ calc_force(MDArena *arena)
        s_calc_angle_force(arena);
        s_calc_dihedral_force(arena);
        s_calc_improper_force(arena);
+#if PIATOM
        s_calc_pibond_force(arena);
+#endif
        s_calc_nonbonded_force(arena);
        s_calc_auxiliary_force(arena);
 

diff --git a/MolLib/MainView.c b/MolLib/MainView.c
index d62cec0..e05422d 100755 (executable)
--- a/MolLib/MainView.c
+++ b/MolLib/MainView.c
@@ -1528,6 +1528,7 @@ skip:
                free(selectFlags);
 }
 
+#if PIATOM
 static void
 drawPiAtoms(MainView *mview)
 {
@@ -1599,6 +1600,7 @@ drawPiAtoms(MainView *mview)
 //     free(vp);
        free(rp);
 }
+#endif
 
 static void
 drawGraphics(MainView *mview)
@@ -1927,7 +1929,9 @@ MainView_drawModel(MainView *mview)
        
        MainViewCallback_clearLabels(mview);
     drawModel(mview);
+#if PIATOM
        drawPiAtoms(mview);
+#endif
        drawUnitCell(mview);
        drawRotationCenter(mview);
        drawGraphics(mview);

diff --git a/MolLib/MolAction.c b/MolLib/MolAction.c
index 85617fe..3cc7134 100644 (file)
--- a/MolLib/MolAction.c
+++ b/MolLib/MolAction.c
@@ -66,11 +66,13 @@ const char *gMolActionClearBox        = "clearBox";
 const char *gMolActionAddParameters   = "addParameters:iGU";
 const char *gMolActionDeleteParameters = "deleteParameters:iG";
 const char *gMolActionAmendBySymmetry = "amendBySymmetry:G;G";
+#if PIATOM
 const char *gMolActionInsertOnePiAtom = "insertOnePiAtom:iCiID";
 const char *gMolActionReplaceOnePiAtom = "replaceOnePiAtom:iCiID";
 const char *gMolActionRemoveOnePiAtom = "removeOnePiAtom:i";
 const char *gMolActionInsertPiBonds   = "insertPiBonds:GI";
 const char *gMolActionRemovePiBonds   = "removePiBonds:G";
+#endif
 
 /*  A Ruby array to retain objects used in MolActionArg  */
 static VALUE sMolActionArgValues = Qfalse;
@@ -1599,6 +1601,7 @@ s_MolActionDeleteParameters(Molecule *mol, MolAction *action, MolAction **actp)
        return 0;
 }
 
+#if PIATOM
 static int
 s_MolActionInsertPiBonds(Molecule *mol, MolAction *action, MolAction **actp)
 {
@@ -1718,6 +1721,7 @@ s_MolActionRemoveOnePiAtom(Molecule *mol, MolAction *action, MolAction **actp)
        MoleculeInvalidatePiConnectionTable(mol);
        return 0;
 }
+#endif  /*  PIATOM  */
 
 int
 MolActionPerform(Molecule *mol, MolAction *action)
@@ -1901,6 +1905,7 @@ MolActionPerform(Molecule *mol, MolAction *action)
                if ((result = s_MolActionDeleteParameters(mol, action, &act2)) != 0)
                        return result;
                needsRebuildMDArena = 1;
+#if PIATOM 
        } else if (strcmp(action->name, gMolActionInsertOnePiAtom) == 0) {
                if ((result = s_MolActionInsertOnePiAtom(mol, action, &act2, 0)) != 0)
                        return result;
@@ -1921,6 +1926,7 @@ MolActionPerform(Molecule *mol, MolAction *action)
                if ((result = s_MolActionRemovePiBonds(mol, action, &act2)) != 0)
                        return result;
                needsRebuildMDArena = 1;
+#endif  /*  PIATOM  */
        } else if (strcmp(action->name, gMolActionNone) == 0) {
                /*  Do nothing  */
        } else {

diff --git a/MolLib/MolAction.h b/MolLib/MolAction.h
index 2da61cf..80c9643 100644 (file)
--- a/MolLib/MolAction.h
+++ b/MolLib/MolAction.h
@@ -67,12 +67,14 @@ extern const char *gMolActionClearBox;
 extern const char *gMolActionAddParameters;
 extern const char *gMolActionDeleteParameters;
 extern const char *gMolActionAmendBySymmetry;
+#if PIATOM
 extern const char *gMolActionInsertOnePiAtom;
 extern const char *gMolActionReplaceOnePiAtom;
 extern const char *gMolActionRemoveOnePiAtom;
 extern const char *gMolActionInsertPiBonds;
 extern const char *gMolActionRemovePiBonds;
-       
+#endif
+
 /*  Special action signatures to invoke the Ruby script. Used as follows:
  *  MolActionCreateAndPerform(mol, SCRIPT_ACTION("vd"), "rotate", vec, angle);
  *    (Will perform 'mol.rotate(vec, angle)')

diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c
index ef84049..ccee809 100755 (executable)
--- a/MolLib/Molecule.c
+++ b/MolLib/Molecule.c
@@ -220,6 +220,7 @@ AtomConnectDeleteEntry(AtomConnect *ac, Int idx)
        AtomConnectResize(ac, ac->count - 1);
 }
 
+#if PIATOM
 void
 PiAtomDuplicate(PiAtom *pa, const PiAtom *cpa)
 {
@@ -245,6 +246,7 @@ PiAtomClean(PiAtom *pa)
                pa->coeffs = NULL;
        }
 }
+#endif
 
 #pragma mark ====== Accessor types ======
 
@@ -367,6 +369,7 @@ MoleculeInitWithMolecule(Molecule *mp2, Molecule *mp)
                NewArray(&(mp2->syms), &(mp2->nsyms), sizeof(Transform), mp->nsyms);
                memmove(mp2->syms, mp->syms, sizeof(Transform) * mp2->nsyms);
        }
+#if PIATOM
        if (mp->npiatoms > 0) {
                NewArray(&(mp2->piatoms), &(mp2->npiatoms), sizeof(PiAtom), mp->npiatoms);
                for (i = 0; i < mp->npiatoms; i++)
@@ -380,6 +383,7 @@ MoleculeInitWithMolecule(Molecule *mp2, Molecule *mp)
                NewArray(&(mp2->piconnects), &(mp2->npiconnects), sizeof(Int), mp->npiconnects);
                memmove(mp2->piconnects, mp->piconnects, sizeof(Int) * mp->npiconnects);
        }
+#endif
        
 /*     mp2->useFlexibleCell = mp->useFlexibleCell; */
        if (mp->nframe_cells > 0) {
@@ -526,6 +530,7 @@ MoleculeClear(Molecule *mp)
                mp->syms = NULL;
                mp->nsyms = 0;
        }
+#if PIATOM
        if (mp->piatoms != NULL) {
                for (i = 0; i < mp->npiatoms; i++) {
                        PiAtomClean(mp->piatoms + i);
@@ -544,6 +549,7 @@ MoleculeClear(Molecule *mp)
                mp->piconnects = NULL;
                mp->npiconnects = 0;
        }
+#endif
        if (mp->selection != NULL) {
                IntGroupRelease(mp->selection);
                mp->selection = NULL;
@@ -1118,6 +1124,7 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char *errbuf, int errbufsi
                                }
                        }
                        continue;
+#if PIATOM
                } else if (strcmp(buf, "!:pi_atoms") == 0) {
                        PiAtom *pp;
                        while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
@@ -1179,6 +1186,7 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char *errbuf, int errbufsi
                pi_atom_constructs_bad_format:
                        snprintf(errbuf, errbufsize, "line %d: bad format in pi_atom_constructs", lineNumber);
                        goto exit;
+#endif  /*  PIATOM  */
                } else if (strcmp(buf, "!:anisotropic_thermal_parameters") == 0) {
                        i = 0;
                        while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
@@ -3932,6 +3940,7 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char *errbuf, int errbu
                fprintf(fp, "\n");
        }
        
+#if PIATOM
        if (mp->npiatoms > 0) {
                PiAtom *pp;
                fprintf(fp, "!:pi_atoms\n");
@@ -3969,7 +3978,8 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char *errbuf, int errbu
                }
                fprintf(fp, "\n");
        }
-       
+#endif  /*  PIATOM  */
+
        if (n_aniso > 0) {
                fprintf(fp, "!:anisotropic_thermal_parameters\n");
                fprintf(fp, "! b11 b22 b33 b12 b13 b23 [sigma; sb11 sb22 sb33 sb12 sb13 sb23]\n");
@@ -5099,6 +5109,7 @@ MoleculeDeserialize(const char *data, Int length, Int *timep)
                        n = len / sizeof(Transform);
                        NewArray(&mp->syms, &mp->nsyms, sizeof(Transform), n);
                        memmove(mp->syms, ptr, len);
+#if PIATOM
                } else if (strcmp(data, "PIATOM") == 0) {
                        const char *ptr2 = ptr + len;
                        mp->piatoms = NULL;
@@ -5127,6 +5138,7 @@ MoleculeDeserialize(const char *data, Int length, Int *timep)
                        n = len / (sizeof(Int) * 4);
                        NewArray(&mp->pibonds, &mp->npibonds, sizeof(Int) * 4, n);
                        memmove(mp->pibonds, ptr, len);
+#endif  /*  PIATOM  */
                } else if (strcmp(data, "TIME") == 0) {
                        if (timep != NULL)
                                *timep = *((Int *)ptr);
@@ -5372,6 +5384,7 @@ MoleculeSerialize(Molecule *mp, Int *outLength, Int *timep)
                len_all += len;
        }
        
+#if PIATOM
        /*  Pi-atoms  */
        if (mp->npiatoms > 0) {
                /*  Estimate the necessary storage first  */
@@ -5421,7 +5434,8 @@ MoleculeSerialize(Molecule *mp, Int *outLength, Int *timep)
                memmove(p, mp->pibonds, sizeof(Int) * 4 * mp->npibonds);
                len_all += len;
        }
-       
+#endif  /*  PIATOM  */
+
        /*  Parameters  */
        if (mp->par != NULL) {
                int type;
@@ -6476,7 +6490,9 @@ MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indice
                return -2;
 
        /*  Create atoms, with avoiding duplicates  */
+#if PIATOM
        MoleculeInvalidatePiConnectionTable(mp);
+#endif
        n0 = n1 = mp->natoms;
        table = (int *)malloc(sizeof(int) * n0);
        if (table == NULL)
@@ -6852,7 +6868,9 @@ MoleculeCreateAnAtom(Molecule *mp, const Atom *ap, int pos)
                        if (mp->impropers[i] >= pos)
                                mp->impropers[i]++;
                }
+#if PIATOM
                MoleculeInvalidatePiConnectionTable(mp);
+#endif
        }
        mp->nframes = -1;  /*  Should be recalculated later  */
        MoleculeIncrementModifyCount(mp);
@@ -6895,7 +6913,9 @@ MoleculeMerge(Molecule *dst, Molecule *src, IntGroup *where, Int resSeqOffset, I
        act = NULL;
 
        __MoleculeLock(dst);
+#if PIATOM
        MoleculeInvalidatePiConnectionTable(dst);
+#endif
        nsrc = src->natoms;
        ndst = dst->natoms;
        if (resSeqOffset < 0)
@@ -7149,6 +7169,7 @@ MoleculeMerge(Molecule *dst, Molecule *src, IntGroup *where, Int resSeqOffset, I
                }
        }
 
+#if PIATOM
        /*  Renumber the existing pi-atoms  */
        if (dst->npiatoms > 0) {
                for (i = 0; i < dst->npiatoms; i++) {
@@ -7221,7 +7242,8 @@ MoleculeMerge(Molecule *dst, Molecule *src, IntGroup *where, Int resSeqOffset, I
                        }
                }
        }
-       
+#endif  /*  PIATOM  */
+
        MoleculeCleanUpResidueTable(dst);
        
        free(new2old);
@@ -7264,7 +7286,9 @@ sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqO
                return 1;  /*  Bad parameter  */
 
        __MoleculeLock(src);
+#if PIATOM
        MoleculeInvalidatePiConnectionTable(src);
+#endif
        
        if (nactions != NULL)
                *nactions = 0;
@@ -7585,7 +7609,8 @@ sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqO
                        }
                }
        }
-       
+
+#if PIATOM
        /*  Copy the pi-atoms  */
        if (src->npiatoms > 0) {
                PiAtom *pp1, *pp2;
@@ -7726,6 +7751,7 @@ sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqO
                
                free(patoms_old2new);
        }
+#endif  /*  PIATOM  */
        
        /*  Clean up  */
        IntGroupRelease(remain_g);
@@ -7829,7 +7855,7 @@ MoleculeExtract(Molecule *src, Molecule **dstp, IntGroup *where, int dummyFlag)
 }
 
 int
-MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds)
+MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds, IntGroup *where, Int autoGenerate)
 {
        int i, j, n1, n2, n;
        Atom *ap;
@@ -7840,39 +7866,17 @@ MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds)
        if (mp->noModifyTopology)
                return -4;  /*  Prohibited operation  */
 
-       /*  Check the bonds  */
-       bonds_tmp = (Int *)malloc(sizeof(Int) * nbonds * 2);
-       if (bonds_tmp == NULL)
+       /*  Note: Duplicates and validity are not checked (the caller must do that)  */
+
+       __MoleculeLock(mp);
+
+       n1 = mp->nbonds;
+       if (AssignArray(&(mp->bonds), &(mp->nbonds), sizeof(Int) * 2, n1 + nbonds - 1, NULL) == NULL
+               || sInsertElementsToArrayAtPositions(mp->bonds, n1, bonds, nbonds, sizeof(Int) * 2, where) != 0) {
+               __MoleculeUnlock(mp);
                return -4;  /*  Out of memory  */
-       n = 0;
-       for (i = 0; i < nbonds; i++) {
-               n1 = bonds[i * 2];
-               n2 = bonds[i * 2 + 1];
-               if (n1 < 0 || n1 >= mp->natoms || n2 < 0 || n2 >= mp->natoms)
-                       return -1;  /*  Bad bond specification  */
-               if (n1 == n2)
-                       return -5;
-               ap = ATOM_AT_INDEX(mp->atoms, n1);
-               /*  Check duplicates  */
-               cp = AtomConnectData(&ap->connect);
-               for (j = 0; j < ap->connect.count; j++) {
-                       if (cp[j] == n2)
-                               break;
-               }
-               if (j == ap->connect.count) {
-                       bonds_tmp[n * 2] = n1;
-                       bonds_tmp[n * 2 + 1] = n2;
-                       n++;
-               }
-       }
-       if (n == 0) {
-               /*  No bonds to add  */
-               free(bonds_tmp);
-               return 0;
        }
        
-       __MoleculeLock(mp);
-
        /*  Add connects[]  */
        for (i = 0; i < n; i++) {
                n1 = bonds_tmp[i * 2];
@@ -7883,106 +7887,96 @@ MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds)
                AtomConnectInsertEntry(&ap->connect, ap->connect.count, n1);
        }
        
-       /*  Expand the array and insert  */
-       n1 = mp->nbonds;
-/*     if (AssignArray(&(mp->bonds), &(mp->nbonds), sizeof(Int) * 2, mp->nbonds + nb - 1, NULL) == NULL
-       || sInsertElementsToArrayAtPositions(mp->bonds, n1, bonds, nb, sizeof(Int) * 2, where) != 0) */
-       if (AssignArray(&(mp->bonds), &(mp->nbonds), sizeof(Int) * 2, mp->nbonds + n - 1, NULL) == NULL)
-               goto panic;
-       memmove(mp->bonds + n1 * 2, bonds_tmp, sizeof(Int) * 2 * n);
-
        /*  Add angles, dihedrals, impropers  */
-       {
-               Int nangles, ndihedrals, nimpropers;
-               Int *angles, *dihedrals, *impropers;
-               Int k, n3, n4;
+       if (autoGenerate) {
+               Int nangles, ndihedrals;
+               Int *angles, *dihedrals;
+               Int k, n3, n4, c1, c2;
                Int *ip, *cp1, *cp2;
                Int temp[4];
-               Atom *ap1, *ap2;
+               Atom *ap1, *ap2, *ap3;
 
-               angles = dihedrals = impropers = NULL;
-               nangles = ndihedrals = nimpropers = 0;
+               angles = dihedrals = NULL;
+               nangles = ndihedrals = 0;
 
                for (i = 0; i < n; i++) {
+                       AtomConnect *ac1, *ac2;
                        n1 = bonds_tmp[i * 2];
                        n2 = bonds_tmp[i * 2 + 1];
                        ap1 = ATOM_AT_INDEX(mp->atoms, n1);
                        ap2 = ATOM_AT_INDEX(mp->atoms, n2);
-                       cp1 = AtomConnectData(&ap1->connect);
-                       cp2 = AtomConnectData(&ap2->connect);
-                       /*  Angles X-n1-n2  */
-                       for (j = 0; j < ap1->connect.count; j++) {
-                               n3 = cp1[j];
-                               if (n3 == n2)
-                                       continue;
-                               temp[0] = n3;
-                               temp[1] = n1;
-                               temp[2] = n2;
-                               for (k = 0; k < nangles; k++) {
-                                       ip = angles + k * 3;
-                                       if (ip[1] == n1 && ((ip[0] == n3 && ip[2] == n2) || (ip[0] == n2 && ip[2] == n3)))
-                                               break;
-                               }
-                               if (k == nangles) {
-                                       if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
-                                               goto panic;
-                               }
-                               /*  Dihedrals X-n1-n2-X  */
-                               for (k = 0; k < ap2->connect.count; k++) {
-                                       n4 = cp2[k];
-                                       if (n4 == n1 || n4 == n3)
-                                               continue;
-                                       temp[3] = n4;
-                                       if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
-                                               goto panic;
+                       if (ap1->anchor == NULL || ap2->anchor == NULL) {
+                               /*  N1-N2-{XY} or N2-N1-{XY} angles (X: connected atom, Y: constitute atom of pi-anchor)  */
+                               for (j = 0; j < 4; j++) {
+                                       switch (j) {
+                                               case 0: temp[0] = n1; temp[1] = n2; ac1 = &ap2->connect; break;  /* N1-N2-X */
+                                               case 1: if (ap2->anchor == NULL) continue; else ac1 = &ap2->anchor->connect; break; /* N1-N2-Y */
+                                               case 2: temp[0] = n2; temp[1] = n1; ac1 = &ap1->connect; break;  /* N2-N1-X */
+                                               case 3: if (ap1->anchor == NULL) continue; else ac1 = &ap1->anchor->connect; break; /* N2-N1-Y */
+                                       }
+                                       cp1 = AtomConnectData(ac1);
+                                       for (k = 0; k < ac1->count; k++) {
+                                               temp[2] = cp1[k];
+                                               if (temp[2] == temp[0])
+                                                       continue;
+                                               if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
+                                                       goto panic;
+                                               /*  Dihedrals N1-N2-X-X or N2-N1-X-X  */
+                                               if (j == 1 || j == 3)
+                                                       continue;
+                                               ap3 = ATOM_AT_INDEX(mp->atoms, temp[2]);
+                                               cp2 = AtomConnectData(&ap3->connect);
+                                               for (kk = 0; kk < ap3->connect.count; kk++) {
+                                                       temp[3] = cp2[kk];
+                                                       if (temp[3] == temp[0] || temp[3] == temp[1])
+                                                               continue;
+                                                       if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
+                                                               goto panic;
+                                               }
+                                       }
                                }
-                               /*  Impropers X-n2-n1-X  */
-                       /*      temp[1] = n2;
-                               temp[2] = n1;
-                               for (k = 0; k < ap1->connect.count; k++) {
-                                       n4 = ap1->connects[k];
-                                       if (n4 == n2 || n4 <= n3)
-                                               continue;
-                                       temp[3] = n4;
-                                       if (AssignArray(&impropers, &nimpropers, sizeof(Int) * 4, nimpropers, temp) == NULL)
-                                               goto panic;
-                               } */
                        }
-                       /*  Angles X-n2-n1  */
-                       for (j = 0; j < ap2->connect.count; j++) {
-                               n3 = cp2[j];
-                               if (n3 == n1)
+                       /*  X-N1-N2-X angles  */
+                       if (ap1->anchor == NULL)
+                               ac1 = &ap1->connect;
+                       else ac1 = &ap1->anchor->connect;
+                       if (ap2->anchor == NULL)
+                               ac2 = &ap2->connect;
+                       else ac2 = &ap2->anchor->connect;
+                       temp[1] = n1;
+                       temp[2] = n2;
+                       cp1 = AtomConnectData(ac1);
+                       cp2 = AtomConnectData(ac2);
+                       for (j = 0; j < ac1->count; j++) {
+                               temp[0] = cp1[j];
+                               if (temp[0] == temp[2])
                                        continue;
-                               temp[0] = n1;
-                               temp[1] = n2;
-                               temp[2] = n3;
-                               for (k = 0; k < nangles; k++) {
-                                       ip = angles + k * 3;
-                                       if (ip[1] == n2 && ((ip[0] == n3 && ip[2] == n1) || (ip[0] == n1 && ip[2] == n3)))
-                                               break;
-                               }
-                               if (k == nangles) {
-                                       if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
+                               if (ATOM_AT_INDEX(mp->atoms, temp[0])->anchor != NULL)
+                                       continue;
+                               for (k = 0; k < ac2->count; k++) {
+                                       temp[3] = cp2[k];
+                                       if (temp[3] == temp[0] || temp[3] == temp[1])
+                                               continue;
+                                       if (ATOM_AT_INDEX(mp->atoms, temp[3])->anchor != NULL)
+                                               continue;
+                                       if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
                                                goto panic;
                                }
                        }
                }
-               temp[0] = kInvalidIndex;
-               if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
-                       goto panic;
-               if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
-                       goto panic;
-               if (AssignArray(&impropers, &nimpropers, sizeof(Int) * 4, nimpropers, temp) == NULL)
-                       goto panic;
-               MoleculeAddAngles(mp, angles, NULL);
-               MoleculeAddDihedrals(mp, dihedrals, NULL);
-               MoleculeAddImpropers(mp, impropers, NULL);
-               if (angles != NULL)
+               temp[0] = kInvalidIndex;  /*  For termination  */
+               if (angles != NULL) {
+                       if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
+                               goto panic;
+                       MoleculeAddAngles(mp, angles, NULL);
                        free(angles);
-               if (dihedrals != NULL)
+               }
+               if (dihedrals != NULL) {
+                       if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
+                               goto panic;
+                       MoleculeAddDihedrals(mp, dihedrals, NULL);
                        free(dihedrals);
-               if (impropers != NULL)
-                       free(impropers);
+               }
        }
        
        MoleculeIncrementModifyCount(mp);
@@ -7996,10 +7990,11 @@ MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds)
        __MoleculeUnlock(mp);
        Panic("Low memory while adding bonds");
        return -1;  /*  Not reached  */
+#endif
 }
 
 int
-MoleculeDeleteBonds(Molecule *mp, Int nbonds, const Int *bonds)
+MoleculeDeleteBonds(Molecule *mp, IntGroup *where, Int **outAutoRemoval, IntGroup **outWhere)
 {
        Int i, j, n1, n2, *cp;
        Atom *ap;
@@ -8482,10 +8477,16 @@ MoleculeFindMissingAngles(Molecule *mol, Int **outAngles)
        angles = NULL;
        for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
                Int *cp = AtomConnectData(&ap->connect);
+               if (ap->anchor != NULL)
+                       continue;
                for (j = 0; j < ap->connect.count; j++) {
                        Int j0 = cp[j];
+                       if (ATOM_AT_INDEX(mol->atoms, j0)->anchor != NULL)
+                               continue;
                        for (k = j + 1; k < ap->connect.count; k++) {
                                Int k0 = cp[k];
+                               if (ATOM_AT_INDEX(mol->atoms, k0)->anchor != NULL)
+                                       continue;
                                if (MoleculeLookupAngle(mol, j0, i, k0) < 0) {
                                        ip = (Int *)AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, NULL);
                                        ip[0] = j0;
@@ -8519,21 +8520,29 @@ MoleculeFindMissingDihedrals(Molecule *mol, Int **outDihedrals)
        dihedrals = NULL;
        for (n2 = 0, ap2 = mol->atoms; n2 < mol->natoms; n2++, ap2 = ATOM_NEXT(ap2)) {
                Int i1, i3, i4, *ip;
+               if (ap2->anchor != NULL)
+                       continue;
                cp2 = AtomConnectData(&ap2->connect);
                for (i3 = 0; i3 < ap2->connect.count; i3++) {
                        n3 = cp2[i3];
                        if (n2 > n3)
                                continue;
                        ap3 = ATOM_AT_INDEX(mol->atoms, n3);
+                       if (ap3->anchor != NULL)
+                               continue;
                        cp3 = AtomConnectData(&ap3->connect);
                        for (i1 = 0; i1 < ap2->connect.count; i1++) {
                                n1 = cp2[i1];
                                if (n1 == n3)
                                        continue;
+                               if (ATOM_AT_INDEX(mol->atoms, n1)->anchor != NULL)
+                                       continue;
                                for (i4 = 0; i4 < ap3->connect.count; i4++) {
                                        n4 = cp3[i4];
                                        if (n2 == n4 || n1 == n4)
                                                continue;
+                                       if (ATOM_AT_INDEX(mol->atoms, n4)->anchor != NULL)
+                                               continue;
                                        if (MoleculeLookupDihedral(mol, n1, n2, n3, n4) < 0) {
                                                ip = (Int *)AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, NULL);
                                                ip[0] = n1;
@@ -8954,8 +8963,9 @@ sMoleculeReorder(Molecule *mp)
        free(old2new);
        free(apArray);
        
+#if PIATOM
        MoleculeInvalidatePiConnectionTable(mp);
-
+#endif
 }
 
 /*  Renumber atoms  */
@@ -9035,6 +9045,7 @@ MoleculeRenumberAtoms(Molecule *mp, const Int *new2old, Int *old2new_out, Int is
                }
        }
        
+#if PIATOM
        if (mp->npiatoms) {
                /*  Renumber the pi-atoms  */
                for (i = 0; i < mp->npiatoms; i++) {
@@ -9053,14 +9064,16 @@ MoleculeRenumberAtoms(Molecule *mp, const Int *new2old, Int *old2new_out, Int is
                                mp->pibonds[i] = old2new[j];
                }
        }
-       
+#endif  /*  PIATOM  */
+
        /*  Renumber the atoms  */
        for (i = 0; i < mp->natoms; i++)
                memmove(ATOM_AT_INDEX(mp->atoms, old2new[i]), ATOM_AT_INDEX(saveAtoms, i), gSizeOfAtomRecord);
        retval = 0;
        
+#if PIATOM
        MoleculeInvalidatePiConnectionTable(mp);
-
+#endif
        MoleculeIncrementModifyCount(mp);
        mp->needsMDRebuild = 1;
 
@@ -9734,6 +9747,7 @@ sMoleculeFragmentSub(Molecule *mp, int idx, IntGroup *result, IntGroup *exatoms)
                        continue;
                sMoleculeFragmentSub(mp, idx2, result, exatoms);
        }
+#if PIATOM
        if (mp->piconnects != NULL) {
                for (i = mp->piconnects[idx]; i < mp->piconnects[idx + 1]; i++) {
                        idx2 = mp->piconnects[i];
@@ -9742,6 +9756,7 @@ sMoleculeFragmentSub(Molecule *mp, int idx, IntGroup *result, IntGroup *exatoms)
                        sMoleculeFragmentSub(mp, idx2, result, exatoms);
                }
        }
+#endif
 }
 
 /*  The molecular fragment (= interconnected atoms) containing the atom n1 and
@@ -9753,7 +9768,9 @@ MoleculeFragmentExcludingAtomGroup(Molecule *mp, int n1, IntGroup *exatoms)
        if (mp == NULL || mp->natoms == 0 || n1 < 0 || n1 >= mp->natoms)
                return NULL;
        result = IntGroupNew();
+#if PIATOM
        MoleculeValidatePiConnectionTable(mp);
+#endif
        sMoleculeFragmentSub(mp, n1, result, exatoms);
        return result;
 }
@@ -9771,7 +9788,9 @@ MoleculeFragmentExcludingAtoms(Molecule *mp, int n1, int argc, int *argv)
        for (i = 0; i < argc; i++)
                IntGroupAdd(exatoms, argv[i], 1);
        result = IntGroupNew();
+#if PIATOM
        MoleculeValidatePiConnectionTable(mp);
+#endif
        sMoleculeFragmentSub(mp, n1, result, exatoms);
        IntGroupRelease(exatoms);
        return result;
@@ -9789,7 +9808,9 @@ MoleculeFragmentWithAtomGroups(Molecule *mp, IntGroup *inatoms, IntGroup *exatom
                return NULL;
        IntGroupIteratorInit(inatoms, &iter);
        result = IntGroupNew();
+#if PIATOM
        MoleculeValidatePiConnectionTable(mp);
+#endif
        while ((i = IntGroupIteratorNext(&iter)) >= 0) {
                sMoleculeFragmentSub(mp, i, result, exatoms);
        }
@@ -9809,7 +9830,9 @@ MoleculeIsFragmentDetachable(Molecule *mp, IntGroup *group, int *n1, int *n2)
        Atom *ap;
        if (mp == NULL || mp->natoms == 0 || group == NULL)
                return 0;  /*  Invalid arguments  */
+#if PIATOM
        MoleculeValidatePiConnectionTable(mp);
+#endif
        bond_count = 0;
        for (i = 0; (i1 = IntGroupGetStartPoint(group, i)) >= 0; i++) {
                i2 = IntGroupGetEndPoint(group, i);
@@ -9827,6 +9850,7 @@ MoleculeIsFragmentDetachable(Molecule *mp, IntGroup *group, int *n1, int *n2)
                                                return 0;  /*  Too many bonds  */
                                }
                        }
+#if PIATOM
                        if (mp->piconnects != NULL) {
                                for (k = mp->piconnects[j]; k < mp->piconnects[j + 1]; k++) {
                                        k2 = mp->piconnects[k];
@@ -9839,6 +9863,7 @@ MoleculeIsFragmentDetachable(Molecule *mp, IntGroup *group, int *n1, int *n2)
                                        }
                                }
                        }
+#endif
                }
        }
        if (bond_count == 1) {
@@ -10260,6 +10285,7 @@ MoleculeFlushFrames(Molecule *mp)
 
 #pragma mark ====== Pi Atoms ======
 
+#if PIATOM
 int
 MoleculeCalculatePiAtomPosition(Molecule *mol, int idx)
 {
@@ -10365,6 +10391,7 @@ MoleculeInvalidatePiConnectionTable(Molecule *mol)
        }
        mol->npiconnects = 0;
 }
+#endif  /*  PIATOM  */
 
 #pragma mark ====== MO calculation ======
 

diff --git a/MolLib/Molecule.h b/MolLib/Molecule.h
index d029a22..d56e706 100755 (executable)
--- a/MolLib/Molecule.h
+++ b/MolLib/Molecule.h
@@ -73,6 +73,12 @@ typedef struct AtomConnect {
        } u;
 } AtomConnect;
 
+typedef struct PiAnchor {
+       AtomConnect connect;
+       Int ncoeffs;
+       Double *coeffs;
+} PiAnchor;
+
 /*  Atom record  */
 typedef struct Atom {
        Int    segSeq;
@@ -100,6 +106,7 @@ typedef struct Atom {
        Vector *frames;
        Symop  symop;    /*  For symmetry-expanded atom  */
        Int    symbase;  /*  The index of original atom for symmetry-expansion  */
+       PiAnchor *anchor;  /*  Non-NULL if this atom is a pi-anchor  */
        Int    labelid;  /*  The label ID; 0 for no label  */
        short  wrap_dx, wrap_dy, wrap_dz; /*  Calculated by md_wrap_coordinates; used only in wrapped output.  */
        Double fix_force; /*  0: no fix, >0: fix at fix_pos with harmonic potential, <0: fix at fix_pos without force  */
@@ -172,7 +179,7 @@ typedef struct XtalCell {
        Double  cellsigma[6];  /*  For crystallographic data; sigma for the cell parameters  */
 } XtalCell;
 
-
+#if PIATOM
 /*  Dummy atoms to represent metal-pi bonds  */
 typedef struct PiAtom {
        char aname[4];
@@ -182,6 +189,7 @@ typedef struct PiAtom {
        Double *coeffs;  /*  The piatom position is given by sum(i, atoms[connect.data[i]] * coeffs[i]) */
        Vector r;        /*  Current position (cache)  */
 } PiAtom;
+#endif
 
 /*  3-Dimensional distribution  */
 typedef struct Cube {
@@ -268,6 +276,8 @@ typedef struct Molecule {
        XtalCell   *cell;
        Int    nsyms;        /*  Symmetry operations; syms are always described in crystallographic units (even when the unit cell is not defined)  */
        Transform *syms;
+
+#if PIATOM
        Int    npiatoms;     /*  Number of "dummy" atoms to represent pi-metal bonds  */
        PiAtom *piatoms;
        Int    npibonds;
@@ -283,7 +293,8 @@ typedef struct Molecule {
                                                     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]]  */
-
+#endif
+       
        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  */
@@ -520,9 +531,11 @@ int MoleculeRemoveFrames(Molecule *mp, IntGroup *group, Vector *outFrame, Vector
 int MoleculeSelectFrame(Molecule *mp, int frame, int copyback);
 int MoleculeFlushFrames(Molecule *mp);
 
+#if PIATOM
 int MoleculeCalculatePiAtomPosition(Molecule *mol, int idx);
 int MoleculeValidatePiConnectionTable(Molecule *mol);
 void MoleculeInvalidatePiConnectionTable(Molecule *mol);
+#endif
        
 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);

diff --git a/MolLib/Ruby_bind/ruby_bind.c b/MolLib/Ruby_bind/ruby_bind.c
index 4ece5e0..4d7c346 100644 (file)
--- a/MolLib/Ruby_bind/ruby_bind.c
+++ b/MolLib/Ruby_bind/ruby_bind.c
@@ -4229,6 +4229,7 @@ s_Molecule_Alloc(VALUE klass)
        return val;
 }
 
+#if PIATOM
 static int
 s_Molecule_AtomOrPiAtomIndexFromValue(Molecule *mol, VALUE val)
 {
@@ -4271,6 +4272,32 @@ s_Molecule_AtomIndexFromValue(Molecule *mol, VALUE val)
        }
        return n;
 }
+#else
+static int
+s_Molecule_AtomIndexFromValue(Molecule *mol, VALUE val)
+{
+       int n, i;
+       char *p;
+       if (FIXNUM_P(val)) {
+               n = FIX2INT(val);
+               if (n >= 0 && n < mol->natoms)
+                       return n;
+               n = -1; /*  No such atom  */
+               val = rb_inspect(val);
+       } else {
+               n = MoleculeAtomIndexFromString(mol, StringValuePtr(val));
+       }
+       if (n >= 0 && n < mol->natoms)
+               return n;
+       if (n == -1)
+               rb_raise(rb_eMolbyError, "no such atom: %s", p);
+       else if (n == -2)
+               rb_raise(rb_eMolbyError, "bad format of atom specification: %s", p);
+       else
+               rb_raise(rb_eMolbyError, "error in converting value to atom index: %s", p);
+       return 0; /* Not reached */
+}
+#endif  /*  _not_ PIATOM  */
 
 static IntGroup *
 s_Molecule_AtomGroupFromValue(VALUE self, VALUE val)
@@ -9428,6 +9455,7 @@ s_Molecule_SearchEquivalentAtoms(int argc, VALUE *argv, VALUE self)
        return val;
 }
 
+#if PIATOM
 /*
  *  call-seq:
  *     create_pi_anchor(name, type, group [, weights]) -> index
@@ -9726,6 +9754,7 @@ s_Molecule_PiAnchorConstructAtIndex(VALUE self, VALUE ival)
        }
        return rb_ary_new4(i, vals);
 }
+#endif  /*  PIATOM  */
 
 /*
  *  call-seq:
@@ -10061,6 +10090,7 @@ Init_Molby(void)
        rb_define_method(rb_cMolecule, "allocate_basis_set_record", s_Molecule_AllocateBasisSetRecord, 3);
        rb_define_method(rb_cMolecule, "search_equivalent_atoms", s_Molecule_SearchEquivalentAtoms, -1);
        
+#if PIATOM
        rb_define_method(rb_cMolecule, "create_pi_anchor", s_Molecule_CreatePiAnchor, -1);
        rb_define_method(rb_cMolecule, "replace_pi_anchor", s_Molecule_ReplacePiAnchor, -1);
        rb_define_method(rb_cMolecule, "insert_pi_anchor", s_Molecule_InsertPiAnchor, -1);
@@ -10074,7 +10104,8 @@ Init_Molby(void)
        rb_define_alias(rb_cMolecule, "remove_pi_anchor_construct", "remove_pi_anchor_constructs");
        rb_define_method(rb_cMolecule, "count_pi_anchor_constructs", s_Molecule_CountPiAnchorConstructs, 0);
        rb_define_method(rb_cMolecule, "pi_anchor_construct", s_Molecule_PiAnchorConstructAtIndex, 1);
-       
+#endif
+
        rb_define_singleton_method(rb_cMolecule, "current", s_Molecule_Current, 0);
        rb_define_singleton_method(rb_cMolecule, "[]", s_Molecule_MoleculeAtIndex, -1);
        rb_define_singleton_method(rb_cMolecule, "open", s_Molecule_Open, -1);

diff --git a/Version b/Version
index 2afc1b4..d4848e0 100644 (file)
--- a/Version
+++ b/Version
@@ -1,2 +1,2 @@
 version = "0.6.4"
-date = "20121025"
+date = "20121026"