OSDN Git Service

Remove previous pi-anchor codes.
authortoshinagata1964 <toshinagata1964@a2be9bc6-48de-4e38-9406-05402d4bc13c>
Mon, 29 Oct 2012 23:43:20 +0000 (23:43 +0000)
committertoshinagata1964 <toshinagata1964@a2be9bc6-48de-4e38-9406-05402d4bc13c>
Mon, 29 Oct 2012 23:43:20 +0000 (23:43 +0000)
git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@312 a2be9bc6-48de-4e38-9406-05402d4bc13c

MolLib/MD/MDCore.c
MolLib/MD/MDForce.c
MolLib/MainView.c
MolLib/MolAction.c
MolLib/MolAction.h
MolLib/Molecule.c
MolLib/Molecule.h
MolLib/Ruby_bind/ruby_bind.c

index 4e22fc8..930afb3 100644 (file)
@@ -293,21 +293,6 @@ 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];
-                       for (ip = results; *ip >= 0; ip++) {
-                               if (n == *ip)
-                                       break;
-                       }
-                       if (*ip < 0) {
-                               *ip++ = n;
-                               *ip = -1;
-                       }
-               }
-       }
-#endif
        for (ip = results; *ip >= 0; ip++)
                ;
        return ip - results;
@@ -562,42 +547,6 @@ 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)
-{
-       Molecule *mol = arena->mol;
-       PiAtom *pp;
-       Int i, j, n1, n2, *cp;
-       for (i = 0; i < mol->npibonds; i++) {
-               if (mol->pibonds[i * 4 + 2] >= 0)
-                       continue;  /*  Not a bond  */
-               n1 = mol->pibonds[i * 4];
-               n2 = mol->pibonds[i * 4 + 1];
-               if (n2 == i1) {
-                       SWAP_INT(n1, n2);
-               } else if (n2 == i2) {
-                       SWAP_INT(n1, n2);
-                       SWAP_INT(i1, i2);
-               } else if (n1 == i2) {
-                       SWAP_INT(i1, i2);
-               } else if (n1 != i1)
-                       continue;  /*  No match with this bond  */
-               /*  Now n1 == i1, so check whether pi-atom "n2" includes i2  */
-               if (n2 < ATOMS_MAX_NUMBER)
-                       continue;
-               pp = mol->piatoms + (n2 - ATOMS_MAX_NUMBER);
-               cp = AtomConnectData(&pp->connect);
-               for (j = pp->connect.count - 1; j >= 0; j--) {
-                       if (cp[j] == i2)
-                               return 1;
-               }
-       }
-       return 0;
-}
-#endif
-
 /*  Find vdw parameters and build the in-use list */
 static int
 s_find_vdw_parameters(MDArena *arena)
@@ -836,13 +785,6 @@ 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;
@@ -930,13 +872,6 @@ 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;
@@ -1025,13 +960,6 @@ 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;
@@ -1198,114 +1126,6 @@ s_find_improper_parameters(MDArena *arena)
        return nmissing;
 }
 
-#if PIATOM
-/*  Find pi-bond parameters  */
-static int
-s_find_pibond_parameters(MDArena *arena)
-{
-       Int i, j, ii, types[4], idx[4], ptype, psize, nmissing = 0;
-       Molecule *mol = arena->mol;
-       UnionPar upar, *up;
-       
-       if (mol->npibonds == 0)
-               return 0;
-       
-       for (i = 0; i < mol->npibonds; i++) {
-               ptype = kDihedralParType;
-               memset(&upar, 0, sizeof(UnionPar));
-               for (j = 0; j < 4; j++) {
-                       ii = mol->pibonds[i * 4 + j];
-                       if (ii < 0) {
-                               if (j < 2)
-                                       ptype = -1;
-                               else if (j == 2)
-                                       ptype = kBondParType;
-                               else
-                                       ptype = kAngleParType;
-                               break;
-                       }
-                       if (ii >= ATOMS_MAX_NUMBER) {
-                               ii -= ATOMS_MAX_NUMBER;
-                               idx[j] = -1;
-                               if (ii >= mol->npiatoms) {
-                                       types[j] = -1;
-                                       ptype = -1;
-                                       break;
-                               } else
-                                       types[j] = mol->piatoms[ii].type;
-                       } else if (ii < mol->natoms) {
-                               idx[j] = ii;
-                               types[j] = ATOM_AT_INDEX(mol->atoms, ii)->type;
-                       } else {
-                               ptype = -1;
-                               break;
-                       }
-               }
-               if (ptype == kBondParType) {
-                       up = (UnionPar *)ParameterLookupBondPar(mol->par, types[0], types[1], idx[0], idx[1], 0);
-                       if (up == NULL)
-                               up = (UnionPar *)ParameterLookupBondPar(gBuiltinParameters, types[0], types[1], idx[0], idx[1], 0);
-                       psize = sizeof(BondPar);
-               } else if (ptype == kAngleParType) {
-                       up = (UnionPar *)ParameterLookupAnglePar(mol->par, types[0], types[1], types[2], idx[0], idx[1], idx[2], 0);
-                       if (up == NULL)
-                               up = (UnionPar *)ParameterLookupAnglePar(gBuiltinParameters, types[0], types[1], types[2], idx[0], idx[1], idx[2], 0);
-                       psize = sizeof(AnglePar);
-               } else if (ptype == kDihedralParType) {
-                       up = (UnionPar *)ParameterLookupDihedralPar(mol->par, types[0], types[1], types[2], types[3], idx[0], idx[1], idx[2], idx[3], 0);
-                       if (up == NULL)
-                               up = (UnionPar *)ParameterLookupDihedralPar(gBuiltinParameters, types[0], types[1], types[2], types[3], idx[0], idx[1], idx[2], idx[3], 0);
-                       psize = sizeof(TorsionPar);
-               } else {
-                       up = NULL;
-                       psize = 0;
-               }
-               if (up != NULL)
-                       memmove(&upar, up, psize);
-               else {
-                       upar.bond.src = -1;
-                       nmissing++;
-               }
-
-               if (arena->debug_result != NULL && arena->debug_output_level > 0) {
-                       char s1[8];
-                       if (i == 0) {
-                               fprintf(arena->debug_result, "\n  Pi anchor parameters\n");
-                               fprintf(arena->debug_result, "  No. atom1 atom2 atom3 atom4 type1  type2  type3  type4  mult r0/a0/phi0 k   per\n");
-                       }
-                       fprintf(arena->debug_result, "%5d ", i + 1);
-                       for (j = 0; j < 4; j++) {
-                               ii = mol->pibonds[i * 4 + j];
-                               if (ii >= 0) {
-                                       fprintf(arena->debug_result, "%5d ", (ii >= ATOMS_MAX_NUMBER ? -(ii - ATOMS_MAX_NUMBER) - 1 : ii));
-                               } else {
-                                       fprintf(arena->debug_result, "----- ");
-                               }
-                       }
-                       for (j = 0; j < 4; j++) {
-                               if (types[j] != -1)
-                                       AtomTypeDecodeToString(types[j], s1);
-                               else
-                                       strcpy(s1, "-----");
-                               fprintf(arena->debug_result, "%-6s ", s1);
-                       }
-                       if (ptype == kBondParType)
-                               fprintf(arena->debug_result, "     %7.3f %7.3f", upar.bond.r0, upar.bond.k*INTERNAL2KCAL);
-                       else if (ptype == kAngleParType)
-                               fprintf(arena->debug_result, "     %7.3f %7.3f", upar.angle.a0*180/PI, upar.angle.k*INTERNAL2KCAL);
-                       else if (ptype == kDihedralParType)
-                               fprintf(arena->debug_result, "%4d %7.3f %7.3f %1d", upar.torsion.mult, upar.torsion.phi0[0]*180/PI, upar.torsion.k[0]*INTERNAL2KCAL, upar.torsion.period[0]);
-                       fprintf(arena->debug_result, "\n");
-               }
-
-               memmove(arena->pi_pars + i, &upar, sizeof(UnionPar));
-       }
-       
-       arena->nmissing += nmissing;
-       return nmissing;
-}
-#endif
-
 /*  Find one fragment, starting from start_index  */
 static void
 s_find_fragment_sub(MDArena *arena, Int start_index, Int fragment_index)
@@ -1755,23 +1575,12 @@ md_prepare(MDArena *arena, int check_only)
                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++) {
@@ -3764,14 +3573,6 @@ 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));
index d084093..2a2e8ea 100644 (file)
@@ -468,93 +468,6 @@ 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)
-{
-       Atom *ap = arena->mol->atoms;
-       Int npibonds = arena->mol->npibonds;
-       Int *pibonds = arena->mol->pibonds;
-       UnionPar *pars = arena->pi_pars;
-       Double energy, val, valk1;
-       Int idx;
-//     Double *energies = &arena->energies[kAngleIndex];
-//     Vector *forces = &arena->forces[kAngleIndex * arena->mol->natoms];
-       Int i;
-       for (i = 0; i < arena->mol->npiatoms; i++) {
-               MoleculeCalculatePiAtomPosition(arena->mol, i);
-       }
-       for (i = 0; i < npibonds; i++, pibonds += 4, pars++) {
-               Vector r[4];
-               Vector f[4], *forces;
-               PiAtom *pp[4];
-               Int j, k, n;
-               for (j = 0; j < 4; j++) {
-                       n = pibonds[j];
-                       if (n < 0)
-                               break;
-                       else if (n < ATOMS_MAX_NUMBER) {
-                               r[j] = ATOM_AT_INDEX(ap, n)->r;
-                               pp[j] = NULL;
-                       } else {
-                               n -= ATOMS_MAX_NUMBER;
-                               pp[j] = arena->mol->piatoms + n;
-                               r[j] = pp[j]->r;
-                       }
-               }
-               if (j == 2) {  /*  Bond  */
-                       Vector r12;
-                       Double w1, w2, k0, k1;
-                       VecSub(r12, r[1], r[0]);
-                       w1 = VecLength(r12);
-                       w2 = w1 - pars->bond.r0;
-                       k0 = pars->bond.k * w2 * w2;         /*  Energy  */
-                       k1 = -2.0 * pars->bond.k * w2 / w1;  /*  Force / r  */
-                       VecScaleSelf(r12, k1);
-                       energy = k0;
-                       VecScale(f[0], r12, -1);
-                       f[1] = r12;
-                       if (arena->debug_result && arena->debug_output_level > 1) {
-                               fprintf(arena->debug_result, "pi bond force %d-%d: r=%f, r0=%f, k1=%f, {%f %f %f}\n", pibonds[0]+1, pibonds[1]+1, w1, pars->bond.r0, k1, r12.x, r12.y, r12.z);
-                       }
-                       idx = kBondIndex;
-               } else if (j == 3) {  /*  Angle  */
-                       if (s_calc_angle_force_one(&r[0], &r[1], &r[2], &pars->angle, &val, &energy, &valk1, &f[0], &f[1], &f[2]) != 0)
-                               continue;
-                       if (arena->debug_result && arena->debug_output_level > 1) {
-                               fprintf(arena->debug_result, "pi angle force %d-%d-%d: a=%f, a0=%f, k1=%f, {%f %f %f}, {%f %f %f}\n", pibonds[0]+1, pibonds[1]+1, pibonds[2]+1, val*180/PI, pars->angle.a0*180/PI, energy, f[0].x, f[0].y, f[0].z, f[1].x, f[1].y, f[1].z);
-                       }
-                       idx = kAngleIndex;
-               } else if (j == 4) {  /*  Dihedral  */
-                       if (s_calc_dihedral_force_one(&r[0], &r[1], &r[2], &r[3], &pars->torsion, &val, &energy, &valk1, &f[0], &f[1], &f[2], &f[3]) != 0)
-                               continue;
-                       if (arena->debug_result && arena->debug_output_level > 1) {
-                               fprintf(arena->debug_result, "pi dihedral force %d-%d-%d-%d: phi=%f, k1=%f, {%f %f %f}, {%f %f %f}, {%f %f %f}\n", pibonds[0]+1, pibonds[1]+1, pibonds[2]+1, pibonds[3]+1, val*180/PI, energy, f[0].x, f[0].y, f[0].z, f[1].x, f[1].y, f[1].z, f[2].x, f[2].y, f[2].z);
-                       }
-                       idx = kDihedralIndex;
-               } else continue; /*  j == 0 or 1? (Cannot happen)  */
-               arena->energies[idx] += energy;
-               forces = &arena->forces[idx * arena->mol->natoms];
-               for (k = 0; k < j; k++) {
-                       n = pibonds[k];
-                       if (n < 0)
-                               break;  /*  Cannot happen  */
-                       else if (n < ATOMS_MAX_NUMBER) {
-                               VecInc(forces[n], f[k]);
-                       } else {
-                               Int *ip = AtomConnectData(&pp[k]->connect);
-                               Int kk;
-                               for (kk = pp[k]->connect.count - 1; kk >= 0; kk--) {
-                                       Int an = ip[kk];  /*  Index of the ring atom  */
-                                       Double co = pp[k]->coeffs[kk];
-                                       VecScaleInc(forces[an], f[k], co);
-                               }
-                       }
-               }
-       }
-}
-#endif  /*  PIATOM  */
-
 /*  ==================================================================== */
 /*  md_check_verlet_list: only for debugging the verlet list generation  */
 /*  ==================================================================== */
@@ -1396,9 +1309,6 @@ 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);
 
index 034994e..ee27660 100755 (executable)
@@ -1555,80 +1555,6 @@ skip:
                free(selectFlags);
 }
 
-#if PIATOM
-static void
-drawPiAtoms(MainView *mview)
-{
-       Int i, j, *cp, nrp;
-       Vector cen;
-       PiAtom *pp;
-       Double rad;
-       GLfloat fval[12];
-       Vector r, *rp;
-       Double d;
-//     static GLfloat sLightGreenColor[] = {0, 1, 0.50, 1};
-       static GLfloat sLightGreenTransColor[] = {0, 1, 0.50, 0.75};
-       Molecule *mol = mview->mol;
-//     Vector *vp = (Vector *)malloc(sizeof(Vector) * mol->npiatoms);
-       nrp = 0;
-       rp = NULL;
-       glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, sLightGreenTransColor);
-       for (i = 0, pp = mol->piatoms; i < mol->npiatoms; i++, pp++) {
-               /*  The pi atom position should be explicitly calculated
-                   (because the base atom(s) may be dragged)  */
-               VecZero(cen);
-               cp = AtomConnectData(&pp->connect);
-               AssignArray(&rp, &nrp, sizeof(Vector), pp->connect.count - 1, NULL);
-               for (j = 0; j < pp->connect.count; j++) {
-                       r = ATOM_AT_INDEX(mol->atoms, cp[j])->r;
-                       if (mview->draggingMode == kMainViewDraggingSelectedAtoms) {
-                               if (MoleculeIsAtomSelected(mol, cp[j]))
-                                       VecInc(r, mview->dragOffset);
-                       }
-                       d = (j < pp->ncoeffs ? pp->coeffs[j] : 0.0);
-                       VecScaleInc(cen, r, d);
-                       rp[j] = r;
-               }
-               pp->r = cen;
-               rad = 0.1;
-               fval[0] = cen.x;
-               fval[1] = cen.y;
-               fval[2] = cen.z;
-               fval[3] = fval[7] = fval[11] = rad;
-               fval[4] = fval[5] = fval[6] = fval[8] = fval[9] = fval[10] = 0.0;
-               drawEllipsoid(fval, fval + 3, fval + 6, fval + 9, 8);
-               for (j = 0; j < pp->connect.count; j++) {
-                       fval[3] = rp[j].x;
-                       fval[4] = rp[j].y;
-                       fval[5] = rp[j].z;
-                       drawCylinder(fval, fval + 3, 0.05, 6, 0);
-               }
-       }
-       for (i = 0, cp = mol->pibonds; i < mol->npibonds; i++, cp += 4) {
-               if (cp[2] >= 0)
-                       continue;  /*  Angle or dihedral  */
-               for (j = 0; j < 2; j++) {
-                       if (cp[j] >= 0 && cp[j] < mol->natoms) {
-                               r = ATOM_AT_INDEX(mol->atoms, cp[j])->r;
-                               if (mview->draggingMode == kMainViewDraggingSelectedAtoms) {
-                                       if (MoleculeIsAtomSelected(mol, cp[j]))
-                                               VecInc(r, mview->dragOffset);
-                               }
-                       } else if (cp[j] >= ATOMS_MAX_NUMBER && cp[j] < ATOMS_MAX_NUMBER + mol->npiatoms) {
-                               r = mol->piatoms[cp[j] - ATOMS_MAX_NUMBER].r;
-                       } else break;
-                       fval[j * 3] = r.x;
-                       fval[j * 3 + 1] = r.y;
-                       fval[j * 3 + 2] = r.z;
-               }
-               if (j == 2)
-                       drawCylinder(fval, fval + 3, 0.05, 6, 0);
-       }
-//     free(vp);
-       free(rp);
-}
-#endif
-
 static void
 drawGraphics(MainView *mview)
 {
@@ -1956,9 +1882,6 @@ MainView_drawModel(MainView *mview)
        
        MainViewCallback_clearLabels(mview);
     drawModel(mview);
-#if PIATOM
-       drawPiAtoms(mview);
-#endif
        drawUnitCell(mview);
        drawRotationCenter(mview);
        drawGraphics(mview);
index f7a83cf..d3d9d6d 100644 (file)
@@ -67,13 +67,6 @@ 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;
@@ -1697,128 +1690,6 @@ s_MolActionDeleteParameters(Molecule *mol, MolAction *action, MolAction **actp)
        return 0;
 }
 
-#if PIATOM
-static int
-s_MolActionInsertPiBonds(Molecule *mol, MolAction *action, MolAction **actp)
-{
-       IntGroup *ig;
-       Int i, j, n, *ip1;
-       ig = action->args[0].u.igval;
-       ip1 = action->args[1].u.arval.ptr;
-       if (ig == NULL || (n = IntGroupGetCount(ig)) == 0)
-               return 0;
-       for (i = 0; i < n; i++) {
-               j = IntGroupGetNthPoint(ig, i);
-               InsertArray(&mol->pibonds, &mol->npibonds, sizeof(Int) * 4, j, 1, ip1 + i * 4);
-       }
-       *actp = MolActionNew(gMolActionRemovePiBonds, ig);
-       MoleculeInvalidatePiConnectionTable(mol);
-       return 0;
-}
-
-static int
-s_MolActionRemovePiBonds(Molecule *mol, MolAction *action, MolAction **actp)
-{
-       IntGroup *ig;
-       Int i, j, n, *ip1;
-       ig = action->args[0].u.igval;
-       if (ig == NULL || (n = IntGroupGetCount(ig)) == 0)
-               return 0;
-       ip1 = (Int *)malloc(sizeof(Int) * 4 * n);
-       for (i = n - 1; i >= 0; i--) {
-               j = IntGroupGetNthPoint(ig, i);
-               DeleteArray(&mol->pibonds, &mol->npibonds, sizeof(Int) * 4, j, 1, ip1 + i * 4);
-       }
-       *actp = MolActionNew(gMolActionInsertPiBonds, ig, n * 4, ip1);
-       free(ip1);
-       MoleculeInvalidatePiConnectionTable(mol);
-       return 0;
-}
-
-static int
-s_MolActionInsertOnePiAtom(Molecule *mol, MolAction *action, MolAction **actp, int willReplace)
-{
-       Int idx, atype, nconnects, *connects, i;
-       char *aname;
-       Double *coeffs;
-       PiAtom *pp;
-       idx = action->args[0].u.ival;
-       if (idx < 0 || idx > mol->npiatoms)
-               return 1;
-       if (willReplace) {
-               if (idx == mol->npiatoms)
-                       return 1;
-               pp = mol->piatoms + idx;
-               nconnects = pp->connect.count;
-               connects = AtomConnectData(&pp->connect);
-               *actp = MolActionNew(gMolActionReplaceOnePiAtom, idx, 4, pp->aname, pp->type, nconnects, connects, pp->ncoeffs, pp->coeffs);
-       } else {
-               *actp = MolActionNew(gMolActionRemoveOnePiAtom, idx);
-       }               
-       aname = action->args[1].u.arval.ptr;
-       atype = action->args[2].u.ival;
-       nconnects = action->args[3].u.arval.nitems;
-       connects = action->args[3].u.arval.ptr;
-       coeffs = action->args[4].u.arval.ptr;
-       if (willReplace) {
-               pp = mol->piatoms + idx;
-               PiAtomClean(pp);
-       } else {
-               pp = (PiAtom *)InsertArray(&mol->piatoms, &mol->npiatoms, sizeof(PiAtom), idx, 1, NULL);
-               memset(pp, 0, sizeof(PiAtom));
-       }
-       memmove(pp->aname, aname, 4);
-       pp->type = atype;
-       AtomConnectResize(&pp->connect, nconnects);
-       memmove(AtomConnectData(&pp->connect), connects, sizeof(Int) * nconnects);
-       NewArray(&pp->coeffs, &pp->ncoeffs, sizeof(Double), nconnects);
-       memmove(pp->coeffs, coeffs, sizeof(Double) * nconnects);
-       if (!willReplace) {
-               for (i = 0; i < mol->npibonds * 4; i++) {
-                       if (mol->pibonds[i] >= ATOMS_MAX_NUMBER + idx)
-                               mol->pibonds[i]++;
-               }
-       }
-       MoleculeInvalidatePiConnectionTable(mol);
-       return 0;
-}
-
-static int
-s_MolActionRemoveOnePiAtom(Molecule *mol, MolAction *action, MolAction **actp)
-{
-       Int idx, i, nconnects, *ip;
-       IntGroup *ig = NULL;
-       PiAtom *pp;
-       idx = action->args[0].u.ival;
-       if (idx < 0 || idx >= mol->npiatoms)
-               return 1;
-       /*  Look for the pibonds entry to be removed  */
-       for (i = 0; i < mol->npibonds * 4; i++) {
-               if (mol->pibonds[i] == idx + ATOMS_MAX_NUMBER) {
-                       if (ig == NULL)
-                               ig = IntGroupNew();
-                       IntGroupAdd(ig, i / 4, 1);
-               }
-       }
-       if (ig != NULL) {
-               MolActionCreateAndPerform(mol, gMolActionRemovePiBonds, ig);
-               IntGroupRelease(ig);
-       }
-       pp = mol->piatoms + idx;
-       nconnects = pp->connect.count;
-       ip = AtomConnectData(&pp->connect);
-       *actp = MolActionNew(gMolActionInsertOnePiAtom, idx, 4, pp->aname, pp->type, nconnects, ip, pp->ncoeffs, pp->coeffs);
-       PiAtomClean(pp);
-       DeleteArray(&mol->piatoms, &mol->npiatoms, sizeof(PiAtom), idx, 1, NULL);
-       for (i = 0; i < mol->npibonds * 4; i++) {
-               if (mol->pibonds[i] > ATOMS_MAX_NUMBER + idx)
-                       mol->pibonds[i]--;
-       }       
-       MoleculeInvalidatePiConnectionTable(mol);
-       return 0;
-}
-#endif  /*  PIATOM  */
-
 int
 MolActionPerform(Molecule *mol, MolAction *action)
 {
@@ -2005,28 +1876,6 @@ 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;
-               needsRebuildMDArena = 1;
-       } else if (strcmp(action->name, gMolActionReplaceOnePiAtom) == 0) {
-               if ((result = s_MolActionInsertOnePiAtom(mol, action, &act2, 1)) != 0)
-                       return result;
-               needsRebuildMDArena = 1;
-       } else if (strcmp(action->name, gMolActionRemoveOnePiAtom) == 0) {
-               if ((result = s_MolActionRemoveOnePiAtom(mol, action, &act2)) != 0)
-                       return result;
-               needsRebuildMDArena = 1;
-       } else if (strcmp(action->name, gMolActionInsertPiBonds) == 0) {
-               if ((result = s_MolActionInsertPiBonds(mol, action, &act2)) != 0)
-                       return result;
-               needsRebuildMDArena = 1;
-       } else if (strcmp(action->name, gMolActionRemovePiBonds) == 0) {
-               if ((result = s_MolActionRemovePiBonds(mol, action, &act2)) != 0)
-                       return result;
-               needsRebuildMDArena = 1;
-#endif  /*  PIATOM  */
        } else if (strcmp(action->name, gMolActionNone) == 0) {
                /*  Do nothing  */
        } else {
index cf7f9cc..70499ee 100644 (file)
@@ -68,13 +68,6 @@ 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);
index 118991d..8f6254a 100755 (executable)
@@ -258,34 +258,6 @@ AtomConnectHasEntry(AtomConnect *ac, Int ent)
        return 0;
 }
 
-#if PIATOM
-void
-PiAtomDuplicate(PiAtom *pa, const PiAtom *cpa)
-{
-       memmove(pa, cpa, sizeof(PiAtom));
-       pa->connect.count = 0;
-       AtomConnectResize(&(pa->connect), cpa->connect.count);
-       memmove(AtomConnectData(&(pa->connect)), AtomConnectData((AtomConnect *)&(cpa->connect)), sizeof(Int) * cpa->connect.count);
-       pa->ncoeffs = 0;
-       pa->coeffs = NULL;
-       if (cpa->ncoeffs > 0) {
-               NewArray(&(pa->coeffs), &(pa->ncoeffs), sizeof(Double), cpa->ncoeffs);
-               memmove(pa->coeffs, cpa->coeffs, sizeof(Double) * cpa->ncoeffs);
-       }
-}
-
-void
-PiAtomClean(PiAtom *pa)
-{
-       AtomConnectResize(&(pa->connect), 0);
-       pa->ncoeffs = 0;
-       if (pa->coeffs != NULL) {
-               free(pa->coeffs);
-               pa->coeffs = NULL;
-       }
-}
-#endif
-
 #pragma mark ====== Accessor types ======
 
 MolEnumerable *
@@ -406,23 +378,8 @@ 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++)
-                       PiAtomDuplicate(mp2->piatoms + i, mp->piatoms + i);
-       }
-       if (mp->npibonds > 0) {
-               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);
-       }
-#endif
-       
-/*     mp2->useFlexibleCell = mp->useFlexibleCell; */
+
+       /*      mp2->useFlexibleCell = mp->useFlexibleCell; */
        if (mp->nframe_cells > 0) {
                if (NewArray(&mp2->frame_cells, &mp2->nframe_cells, sizeof(Vector) * 4, mp->nframe_cells) == NULL)
                        goto error;
@@ -566,26 +523,6 @@ 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);
-               }
-               free(mp->piatoms);
-               mp->piatoms = NULL;
-               mp->npiatoms = 0;
-       }
-       if (mp->pibonds != NULL) {
-               free(mp->pibonds);
-               mp->pibonds = NULL;
-               mp->npibonds = 0;
-       }
-       if (mp->piconnects != NULL) {
-               free(mp->piconnects);
-               mp->piconnects = NULL;
-               mp->npiconnects = 0;
-       }
-#endif
        if (mp->selection != NULL) {
                IntGroupRelease(mp->selection);
                mp->selection = NULL;
@@ -1245,69 +1182,6 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf)
                                }
                        }
                        continue;
-#if PIATOM
-               } else if (strcmp(buf, "!:pi_atoms") == 0) {
-                       PiAtom *pp;
-                       while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
-                               if (buf[0] == '!')
-                                       continue;
-                               if (buf[0] == '\n')
-                                       break;
-                               if (sscanf(buf, "%6s %6s %d", cbuf[0], cbuf[1], &ibuf[0]) < 3) {
-                                       s_append_asprintf(errbuf, "line %d: pi atoms info cannot be read", lineNumber);
-                                       goto err_exit;
-                               }
-                               pp = (PiAtom *)AssignArray(&mp->piatoms, &mp->npiatoms, sizeof(PiAtom), mp->npiatoms, NULL);
-                               memset(pp, 0, sizeof(PiAtom));
-                               strncpy(pp->aname, cbuf[0], 4);
-                               pp->type = AtomTypeEncodeToUInt(cbuf[1]);
-                               if (ibuf[0] <= 0)
-                                       continue;
-                               AtomConnectResize(&pp->connect, ibuf[0]);
-                               ip = AtomConnectData(&pp->connect);
-                               NewArray(&pp->coeffs, &pp->ncoeffs, sizeof(Double), ibuf[0]);
-                               for (i = 0; i < ibuf[0]; i++) {
-                                       if (ReadLine(buf, sizeof buf, fp, &lineNumber) <= 0) {
-                                               s_append_asprintf(errbuf, "line %d: unexpected end of file during reading pi atoms info", lineNumber);
-                                               goto err_exit;
-                                       }
-                                       if (sscanf(buf, "%d %lf", &ibuf[1], &dbuf[0]) < 2) {
-                                               s_append_asprintf(errbuf, "line %d: bad format during pi atoms info", lineNumber);
-                                               goto err_exit;
-                                       }
-                                       ip[i] = ibuf[1];
-                                       pp->coeffs[i] = dbuf[0];
-                               }
-                       }
-                       continue;
-               } else if (strcmp(buf, "!:pi_atom_constructs") == 0) {
-                       while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
-                               if (buf[0] == '!')
-                                       continue;
-                               if (buf[0] == '\n')
-                                       break;
-                               /* a1 b1 c1 d1 a2 b2 c2 d2 */ 
-                               i = sscanf(buf, "%d %d %d %d %d %d %d %d", &ibuf[0], &ibuf[1], &ibuf[2], &ibuf[3], &ibuf[4], &ibuf[5], &ibuf[6], &ibuf[7]);
-                               if (i == 0 || i % 4 != 0)
-                                       goto pi_atom_constructs_bad_format;
-                               for (j = 0; j < i; j++) {
-                                       if (ibuf[j] <= -2) {
-                                               ibuf[j] = (-ibuf[j] - 2) + ATOMS_MAX_NUMBER;
-                                               if (ibuf[j] - ATOMS_MAX_NUMBER >= mp->npiatoms)
-                                                       goto pi_atom_constructs_bad_format;
-                                       } else if (ibuf[j] >= mp->natoms) {
-                                               goto pi_atom_constructs_bad_format;
-                                       }
-                                       if (j % 4 == 3) {
-                                               AssignArray(&mp->pibonds, &mp->npibonds, sizeof(Int) * 4, mp->npibonds, &ibuf[j - 3]);
-                                       }
-                               }
-                       }
-                       continue;
-               pi_atom_constructs_bad_format:
-                       s_append_asprintf(errbuf, "line %d: bad format in pi_atom_constructs", lineNumber);
-                       goto err_exit;
-#endif  /*  PIATOM  */
                } else if (strcmp(buf, "!:anisotropic_thermal_parameters") == 0) {
                        i = 0;
                        while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
@@ -4087,46 +3961,6 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf)
                fprintf(fp, "\n");
        }
        
-#if PIATOM
-       if (mp->npiatoms > 0) {
-               PiAtom *pp;
-               fprintf(fp, "!:pi_atoms\n");
-               fprintf(fp, "! name type n; a1 coeff1; a2 coeff2; ...; a_n coeff_n\n");
-               for (i = 0, pp = mp->piatoms; i < mp->npiatoms; i++, pp++) {
-                       strncpy(bufs[0], pp->aname, 4);
-                       bufs[0][4] = 0;
-                       AtomTypeDecodeToString(pp->type, bufs[1]);
-                       bufs[1][6] = 0;
-                       for (j = 0; j < 2; j++) {
-                               if (bufs[j][0] == 0) {
-                                       bufs[j][0] = '_';
-                                       bufs[j][1] = 0;
-                               }
-                       }
-                       fprintf(fp, "%s %s %d\n", bufs[0], bufs[1], pp->connect.count);
-                       ip = AtomConnectData(&pp->connect);
-                       for (j = 0; j < pp->connect.count; j++) {
-                               fprintf(fp, "%d %g\n", ip[j], (j < pp->ncoeffs ? pp->coeffs[j] : 0.0));
-                       }
-               }
-               fprintf(fp, "\n");
-       }
-       
-       if (mp->npibonds > 0) {
-               fprintf(fp, "!:pi_atom_constructs\n");
-               fprintf(fp, "! a1 b1 c1 d1 a2 b2 c2 d2\n");
-               for (i = 0; i < mp->npibonds * 4; i++) {
-                       j = mp->pibonds[i];
-                       if (j >= ATOMS_MAX_NUMBER)
-                               j = -2 - (j - ATOMS_MAX_NUMBER);
-                       else if (j < 0)
-                               j = -1;
-                       fprintf(fp, "%d%c", j, (i % 8 == 7 || i == mp->npibonds * 4 - 1 ? '\n' : ' '));
-               }
-               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");
@@ -5271,36 +5105,6 @@ MoleculeDeserialize(const char *data, Int length, Int *timep)
                                }
                                ptr += sizeof(Int) * (2 + n) + sizeof(Double) * n;
                        }
-#if PIATOM
-               } else if (strcmp(data, "PIATOM") == 0) {
-                       const char *ptr2 = ptr + len;
-                       mp->piatoms = NULL;
-                       mp->npiatoms = 0;
-                       while (ptr < ptr2) {
-                               PiAtom pa;
-                               memset(&pa, 0, sizeof(pa));
-                               n = offsetof(PiAtom, connect);
-                               memmove(&pa, ptr, n);
-                               ptr += n;
-                               n = *((Int *)ptr);
-                               if (n > 0) {
-                                       AtomConnectResize(&pa.connect, n);
-                                       memmove(AtomConnectData(&pa.connect), ptr + sizeof(Int), n * sizeof(Int));
-                               }
-                               ptr += sizeof(Int) * (n + 1);
-                               n = *((Int *)ptr);
-                               if (n > 0) {
-                                       NewArray(&pa.coeffs, &pa.ncoeffs, sizeof(Double), n);
-                                       memmove(pa.coeffs, ptr + sizeof(Int), n * sizeof(Double));
-                               }
-                               ptr += sizeof(Int) + sizeof(Double) * n;
-                               AssignArray(&mp->piatoms, &mp->npiatoms, sizeof(PiAtom), mp->npiatoms, &pa);
-                       }
-               } else if (strcmp(data, "PIBOND") == 0) {
-                       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);
@@ -5586,58 +5390,6 @@ MoleculeSerialize(Molecule *mp, Int *outLength, Int *timep)
                len_all += len;
        }
        
-#if PIATOM
-       /*  Pi-atoms  */
-       if (mp->npiatoms > 0) {
-               /*  Estimate the necessary storage first  */
-               len = 8 + sizeof(Int);
-               for (i = 0; i < mp->npiatoms; i++) {
-                       len += offsetof(PiAtom, connect);  /*  Members before 'connect' is stored as they are  */
-                       len += sizeof(Int) * (1 + mp->piatoms[i].connect.count); /* Array of Int's */
-                       len += sizeof(Int) + sizeof(Double) * mp->piatoms[i].ncoeffs; /* Array of Double's */
-               }
-               ptr = (char *)realloc(ptr, len_all + len);
-               if (ptr == NULL)
-                       goto out_of_memory;
-               p = ptr + len_all;
-               memmove(p, "PIATOM\0\0", 8);
-               *((Int *)(p + 8)) = len - (8 + sizeof(Int));
-               p += 8 + sizeof(Int);
-               for (i = 0; i < mp->npiatoms; i++) {
-                       int len0;
-                       PiAtom *pp = &(mp->piatoms[i]);
-                       len0 = offsetof(PiAtom, connect);
-                       memmove(p, pp, len0);
-                       p += len0;
-                       len0 = pp->connect.count * sizeof(Int);
-                       *((Int *)p) = pp->connect.count;
-                       if (len0 > 0)
-                               memmove(p + sizeof(Int), AtomConnectData(&(pp->connect)), len0);
-                       p += sizeof(Int) + len0;
-                       len0 = pp->ncoeffs * sizeof(Double);
-                       *((Int *)p) = pp->ncoeffs;
-                       if (len0 > 0)
-                               memmove(p + sizeof(Int), pp->coeffs, len0);
-                       p += sizeof(Int) + len0;
-               }
-               len_all += len;
-       }
-       
-       /*  Pi-atom constructs  */
-       if (mp->npibonds > 0) {
-               len = 8 + sizeof(Int) + sizeof(Int) * 4 * mp->npibonds;
-               ptr = (char *)realloc(ptr, len_all + len);
-               if (ptr == NULL)
-                       goto out_of_memory;
-               p = ptr + len_all;
-               memmove(p, "PIBOND\0\0", 8);
-               *((Int *)(p + 8)) = sizeof(Int) * 4 * mp->npibonds;
-               p += 8 + sizeof(Int);
-               memmove(p, mp->pibonds, sizeof(Int) * 4 * mp->npibonds);
-               len_all += len;
-       }
-#endif  /*  PIATOM  */
-
        /*  Parameters  */
        if (mp->par != NULL) {
                int type;
@@ -6687,9 +6439,6 @@ 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)
@@ -7101,9 +6850,6 @@ 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);
@@ -7229,9 +6975,7 @@ 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)
@@ -7487,81 +7231,6 @@ 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++) {
-                       PiAtom *pp;
-                       pp = &dst->piatoms[i];
-                       cp = AtomConnectData(&pp->connect);
-                       for (j = 0; j < pp->connect.count; j++) {
-                               cp[j] = old2new[cp[j]];
-                       }
-               }
-               if (dst->npibonds > 0) {
-                       cp = dst->pibonds;
-                       for (i = 0; i < dst->npibonds * 4; i++) {
-                               if (cp[i] < 0)
-                                       continue;
-                               else if (cp[i] >= ATOMS_MAX_NUMBER)
-                                       continue;
-                               else {
-                                       cp[i] = old2new[cp[i]];
-                               }
-                       }
-               }
-       }
-       
-       /*  Copy the pi-atoms  */
-       if (src->npiatoms > 0 && forUndo == 0) {
-               int nsrcp, ndstp;
-               nsrcp = src->npiatoms;
-               ndstp = dst->npiatoms;
-               if (AssignArray(&dst->piatoms, &dst->npiatoms, sizeof(PiAtom), nsrcp + ndstp - 1, NULL) == NULL)
-                       goto panic;
-               for (i = 0; i < nsrcp; i++) {
-                       PiAtom *pp;
-                       pp = &dst->piatoms[ndstp + i];
-                       PiAtomDuplicate(pp, &src->piatoms[i]);
-                       cp = AtomConnectData(&pp->connect);
-                       for (j = 0; j < pp->connect.count; j++) {
-                               cp[j] = old2new[ndst + cp[j]];
-                       }
-                       if (nactions != NULL) {
-                               /*  This is very inefficient, yet should not cause big problems
-                                   because the number of piatoms in the molecule is usually limited.  */
-                               act = MolActionNew(gMolActionRemoveOnePiAtom, ndstp + i);
-                               AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
-                               act = NULL;
-                       }
-               }
-               if (src->npibonds > 0) {
-                       n1 = src->npibonds;
-                       n2 = dst->npibonds;
-                       if (AssignArray(&dst->pibonds, &dst->npibonds, sizeof(Int) * 4, n1 + n2 - 1, NULL) == NULL)
-                               goto panic;
-                       cp = &dst->pibonds[n2 * 4];
-                       memmove(cp, src->pibonds, sizeof(Int) * 4 * n1);
-                       for (i = 0; i < 4 * n1; i++) {
-                               /*  Renumber the pi-atom constructs  */
-                               if (cp[i] < 0)
-                                       continue;
-                               else if (cp[i] < ATOMS_MAX_NUMBER)
-                                       cp[i] = old2new[ndst + cp[i]];  /*  Ordinary atoms  */
-                               else
-                                       cp[i] += ndstp;  /*  pi-atoms  */
-                       }
-                       if (nactions != NULL) {
-                               ig = IntGroupNewWithPoints(n2, n1, -1);
-                               act = MolActionNew(gMolActionRemovePiBonds, ig);
-                               AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
-                               act = NULL;
-                               IntGroupRelease(ig);
-                       }
-               }
-       }
-#endif  /*  PIATOM  */
-
        MoleculeCleanUpResidueTable(dst);
        
        free(new2old);
@@ -7606,9 +7275,6 @@ sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqO
                return 1;  /*  Bad parameter  */
 
        __MoleculeLock(src);
-#if PIATOM
-       MoleculeInvalidatePiConnectionTable(src);
-#endif
        
        act = NULL;
        
@@ -8000,149 +7666,6 @@ sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqO
                }
        }
 
-#if PIATOM
-       /*  Copy the pi-atoms  */
-       if (src->npiatoms > 0) {
-               PiAtom *pp1, *pp2;
-               Int *patoms_old2new;
-               n1 = src->npiatoms;
-               patoms_old2new = (Int *)calloc(sizeof(Int), src->npiatoms);
-               if (patoms_old2new == NULL)
-                       goto panic;
-               for (i = 0, pp1 = src->piatoms; i < src->npiatoms; i++, pp1++) {
-                       /*  Is this entry to be copied to dst?  */
-                       cp = AtomConnectData(&pp1->connect);
-                       for (j = pp1->connect.count - 1; j >= 0; j--) {
-                               if (old2new[cp[j]] < nsrcnew)
-                                       break;
-                       }
-                       if (j < 0) {
-                               /*  Copy this entry  */
-                               patoms_old2new[i] = dst->npiatoms;
-                               pp2 = AssignArray(&dst->piatoms, &dst->npiatoms, sizeof(PiAtom), dst->npiatoms, NULL);
-                               PiAtomDuplicate(pp2, pp1);
-                               cp = AtomConnectData(&pp2->connect);
-                               for (j = 0; j < pp2->connect.count; j++)
-                                       cp[j] = old2new[cp[j]] - nsrcnew;
-                       } else {
-                               patoms_old2new[i] = -1;
-                       }
-               }
-               /*  Copy the piatom constructs to dst  */
-               for (i = 0; i < src->npibonds; i++) {
-                       for (j = 0; j < 4; j++) {
-                               n2 = src->pibonds[i * 4 + j];
-                               if (n2 >= ATOMS_MAX_NUMBER) {
-                                       if (patoms_old2new[n2 - ATOMS_MAX_NUMBER] < 0)
-                                               break;
-                               } else if (n2 >= 0) {
-                                       if (old2new[n2] < nsrcnew)
-                                               break;
-                               }
-                       }
-                       if (j >= 4) {
-                               /*  Copy this entry  */
-                               cp = (Int *)AssignArray(&dst->pibonds, &dst->npibonds, sizeof(Int) * 4, dst->npibonds, &src->pibonds[i * 4]);
-                               for (j = 0; j < 4; j++) {
-                                       if (cp[j] >= ATOMS_MAX_NUMBER) {
-                                               cp[j] = patoms_old2new[cp[j] - ATOMS_MAX_NUMBER] + ATOMS_MAX_NUMBER;
-                                       } else if (cp[j] >= 0) {
-                                               cp[j] = old2new[cp[j]] - nsrcnew;
-                                       }
-                               }
-                       }
-               }
-               if (moveFlag) {
-                       Int npibonds_to_go, *pibonds_to_go;
-                       IntGroup *pibonds_group;
-
-                       /*  Remove the piatom entries containing non-remaining atoms. Note: the piatom
-                        entries that do not remain in src and not copied to dst will disappear.  */
-                       n2 = 0;
-                       for (i = 0; i < src->npiatoms; i++) {
-                               /*  Is this entry to be removed?  */
-                               pp1 = &src->piatoms[i];
-                               cp = AtomConnectData(&pp1->connect);
-                               for (j = pp1->connect.count - 1; j >= 0; j--) {
-                                       if (old2new[cp[j]] >= nsrcnew)
-                                               break;
-                               }
-                               patoms_old2new[i] = (j < 0 ? n2++ : -1);
-                       }
-
-                       /*  Remove pibonds first (necessary for undo)  */
-                       npibonds_to_go = 0;
-                       pibonds_to_go = NULL;
-                       pibonds_group = NULL;
-                       for (i = src->npibonds - 1; i >= 0; i--) {
-                               for (j = 0; j < 4; j++) {
-                                       n3 = src->pibonds[i * 4 + j];
-                                       if (n3 >= ATOMS_MAX_NUMBER) {
-                                               if (patoms_old2new[n3 - ATOMS_MAX_NUMBER] < 0)
-                                                       break;
-                                       } else if (n3 >= 0) {
-                                               if (old2new[n3] >= nsrcnew)
-                                                       break;
-                                       }
-                               }
-                               if (j < 4) {
-                                       /*  Remove  */
-                                       if (nactions != NULL) {
-                                               /*  Since we are scanning pibonds[] from the end, the new entry should be inserted
-                                                   at the top of the pibonds_to_go[] array.  */
-                                               InsertArray(&pibonds_to_go, &npibonds_to_go, sizeof(Int) * 4, 0, 1, src->pibonds + i * 4);
-                                               if (pibonds_group == NULL)
-                                                       pibonds_group = IntGroupNew();
-                                               IntGroupAdd(pibonds_group, i, 1);
-                                       }
-                                       DeleteArray(&src->pibonds, &src->npibonds, sizeof(Int) * 4, i, 1, NULL);
-                               } else {
-                                       /*  Renumber  */
-                                       cp = &src->pibonds[i * 4];
-                                       for (j = 0; j < 4; j++) {
-                                               n3 = cp[j];
-                                               if (n3 >= ATOMS_MAX_NUMBER)
-                                                       cp[j] = patoms_old2new[n3 - ATOMS_MAX_NUMBER] + ATOMS_MAX_NUMBER;
-                                               else if (n3 >= 0)
-                                                       cp[j] = old2new[n3];
-                                       }
-                               }
-                       }
-                       if (nactions != NULL && pibonds_to_go != NULL) {
-                               act = MolActionNew(gMolActionInsertPiBonds, pibonds_group, npibonds_to_go * 4, pibonds_to_go);
-                               AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
-                               act = NULL;
-                               free(pibonds_to_go);
-                       }
-                       if (pibonds_group != NULL)
-                               IntGroupRelease(pibonds_group);
-
-                       for (i = src->npiatoms - 1; i >= 0; i--) {
-                               pp1 = &src->piatoms[i];
-                               if (patoms_old2new[i] < 0) {
-                                       /*  Remove the entries  */
-                                       /*  (If forUndo is true, these entries should already have been removed.)  */
-                                       if (nactions != NULL) {
-                                               act = MolActionNew(gMolActionInsertOnePiAtom, i, 4, pp1->aname, pp1->type, pp1->connect.count, AtomConnectData(&pp1->connect), pp1->ncoeffs, pp1->coeffs);
-                                               AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
-                                               act = NULL;
-                                       }
-                                       PiAtomClean(pp1);
-                                       DeleteArray(&src->piatoms, &src->npiatoms, sizeof(PiAtom), i, 1, NULL);
-                               } else {
-                                       /*  Renumber  */
-                                       cp = AtomConnectData(&pp1->connect);
-                                       for (j = 0; j < pp1->connect.count; j++) {
-                                               cp[j] = old2new[cp[j]];
-                                       }
-                               }
-                       }
-               }
-               
-               free(patoms_old2new);
-       }
-#endif  /*  PIATOM  */
-       
        /*  Clean up  */
        IntGroupRelease(remain_g);
        MoleculeCleanUpResidueTable(src);
@@ -9392,10 +8915,6 @@ sMoleculeReorder(Molecule *mp)
        free(newAtoms);
        free(old2new);
        free(apArray);
-       
-#if PIATOM
-       MoleculeInvalidatePiConnectionTable(mp);
-#endif
 }
 
 /*  Renumber atoms  */
@@ -9475,35 +8994,11 @@ 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++) {
-                       PiAtom *pp = &mp->piatoms[i];
-                       Int *cp = AtomConnectData(&pp->connect);
-                       for (j = 0; j < pp->connect.count; j++)
-                               cp[j] = old2new[cp[j]];
-               }
-       }
-       
-       if (mp->npibonds) {
-               /*  Renumber the pi-atom constructs  */
-               for (i = 0; i < mp->npibonds * 4; i++) {
-                       j = mp->pibonds[i];
-                       if (j >= 0 && j < mp->natoms)
-                               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;
 
@@ -10188,16 +9683,6 @@ sMoleculeFragmentSub(Molecule *mp, int idx, IntGroup *result, IntGroup *exatoms)
                        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];
-                       if (IntGroupLookup(result, idx2, NULL))
-                               continue;
-                       sMoleculeFragmentSub(mp, idx2, result, exatoms);
-               }
-       }
-#endif
 }
 
 /*  The molecular fragment (= interconnected atoms) containing the atom n1 and
@@ -10209,9 +9694,6 @@ 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;
 }
@@ -10229,9 +9711,6 @@ 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;
@@ -10249,9 +9728,6 @@ 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);
        }
@@ -10271,9 +9747,6 @@ 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);
@@ -10305,20 +9778,6 @@ MoleculeIsFragmentDetachable(Molecule *mp, IntGroup *group, int *n1, int *n2)
                                        }
                                }                                       
                        }
-#if PIATOM
-                       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  */
-                                       }
-                               }
-                       }
-#endif
                }
        }
        if (bond_count == 1) {
@@ -10740,8 +10199,6 @@ MoleculeFlushFrames(Molecule *mp)
 
 #pragma mark ====== Pi Atoms ======
 
-#if !defined(PIATOM)
-
 void
 MoleculeCalculatePiAnchorPosition(Molecule *mol, int idx)
 {
@@ -10892,116 +10349,6 @@ MoleculeSetPiAnchorList(Molecule *mol, Int idx, Int nentries, Int *entries, Doub
        return 0;
 }
 
-#endif
-
-#if PIATOM
-int
-MoleculeCalculatePiAtomPosition(Molecule *mol, int idx)
-{
-       PiAtom *pp;
-       Int i, *cp;
-       if (mol == NULL || idx < 0 || idx >= mol->npiatoms)
-               return -1;
-       pp = mol->piatoms + idx;
-       cp = AtomConnectData(&pp->connect);
-       VecZero(pp->r);
-       for (i = pp->connect.count - 1; i >= 0; i--) {
-               Vector rr = ATOM_AT_INDEX(mol->atoms, cp[i])->r;
-               VecScaleInc(pp->r, rr, pp->coeffs[i]);
-       }
-       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;
-}
-#endif  /*  PIATOM  */
-
 #pragma mark ====== MO calculation ======
 
 /*  Calculate an MO value for a single point.  */
index 34de218..9e21896 100755 (executable)
@@ -179,18 +179,6 @@ 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];
-       UInt type;
-       AtomConnect connect;
-       Int  ncoeffs;
-       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 {
        Int idn;             /*  Integer identifier (such as MO number)  */
@@ -277,24 +265,6 @@ typedef struct Molecule {
        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;
-       Int    *pibonds;     /*  Array to represent bond/angle/dihedral including piatoms. */
-                         /* [n1, n2, -1, -1]: bonds,
-                                                       [n1, n2, n3, -1]: angle,
-                                                   [n1, n2, n3, n4]: dihedral,
-                                                   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]]  */
-#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  */
@@ -537,16 +507,8 @@ int MoleculeRemoveFrames(Molecule *mp, IntGroup *group, Vector *outFrame, Vector
 int MoleculeSelectFrame(Molecule *mp, int frame, int copyback);
 int MoleculeFlushFrames(Molecule *mp);
 
-#if !defined(PIATOM)
 void MoleculeCalculatePiAnchorPosition(Molecule *mol, int idx);
 int MoleculeSetPiAnchorList(Molecule *mol, Int idx, Int nentries, Int *entries, Double *weights, Int *nUndoActions, struct MolAction ***undoActions);
-#endif
-       
-#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);
index dccf384..68b4043 100644 (file)
@@ -4331,50 +4331,6 @@ s_Molecule_Alloc(VALUE klass)
        return val;
 }
 
-#if PIATOM
-static int
-s_Molecule_AtomOrPiAtomIndexFromValue(Molecule *mol, VALUE val)
-{
-       int n, i;
-       char *p;
-       if (FIXNUM_P(val)) {
-               n = FIX2INT(val);
-               if (n >= 0 && n < mol->natoms)
-                       return n;
-               else if (n < 0 && -n - 1 < mol->npiatoms)
-                       return ATOMS_MAX_NUMBER + (-n - 1);  /*  PiAtom  */
-               n = -1; /*  No such atom  */
-               val = rb_inspect(val);
-       } else {
-               n = MoleculeAtomIndexFromString(mol, StringValuePtr(val));
-       }
-       if (n >= 0 && n < mol->natoms)
-               return n;
-       p = StringValuePtr(val);
-       for (i = 0; i < mol->npiatoms; i++) {
-               if (strcmp(p, mol->piatoms[i].aname) == 0)
-                       return ATOMS_MAX_NUMBER + i;  /*  PiAtom; should be looked up after normal atoms  */
-       }
-       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 */
-}
-
-static int
-s_Molecule_AtomIndexFromValue(Molecule *mol, VALUE val)
-{
-       int n = s_Molecule_AtomOrPiAtomIndexFromValue(mol, val);
-       if (n >= ATOMS_MAX_NUMBER) {
-               val = rb_inspect(val);
-               rb_raise(rb_eMolbyError, "no such atom: %s", StringValuePtr(val));
-       }
-       return n;
-}
-#else
 static int
 s_Molecule_AtomIndexFromValue(Molecule *mol, VALUE val)
 {
@@ -4399,7 +4355,6 @@ s_Molecule_AtomIndexFromValue(Molecule *mol, VALUE val)
                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)
@@ -9541,7 +9496,6 @@ s_Molecule_SearchEquivalentAtoms(int argc, VALUE *argv, VALUE self)
        return val;
 }
 
-#if !defined(PIATOM)
 /*
  *  call-seq:
  *     create_pi_anchor(name, group [, type] [, weights] [, index]) -> AtomRef
@@ -9641,309 +9595,6 @@ s_Molecule_CreatePiAnchor(int argc, VALUE *argv, VALUE self)
     return Data_Wrap_Struct(rb_cAtomRef, 0, (void (*)(void *))AtomRefRelease, aref);
 }
 
-#endif
-
-#if PIATOM
-/*
- *  call-seq:
- *     create_pi_anchor(name, type, group [, weights]) -> index
- *     insert_pi_anchor(index, name, type, group [, weights]) -> index
- *     replace_pi_anchor(index, name, type, group [, weights]) -> index
- *
- *  Create or replace a "pi anchor", which is an anchor point to define a metal-pi bond.
- *  If index is negative or no less than the current number of pi anchors, create
- *  a new anchor. Otherwise, replace the existing anchor.
- *  Name and type are Strings, and are similar to those in atoms.
- *  Group is a group of atoms to define a pi system to be bound to the metal.
- *  Weights (optional) is an Array of Floats, which determine the significance
- *  of the component atoms. If not given, then 1.0/N (N is the number of atoms
- *  in the group) is assumed for all atoms.
- */
-static VALUE
-s_Molecule_CreateOrReplacePiAnchor(int func, int argc, VALUE *argv, VALUE self)
-{
-       Molecule *mol;
-       Int idx, i, j, atype, *ip;
-       VALUE ival, nval, tval, gval, wval;
-       char *np, *tp, aname[6];
-       Double *dp;
-       IntGroup *ig;
-    Data_Get_Struct(self, Molecule, mol);
-       if (func == 0) { /*  create  */
-               rb_scan_args(argc, argv, "31", &nval, &tval, &gval, &wval);
-               idx = mol->npiatoms;
-       } else {
-               rb_scan_args(argc, argv, "41", &ival, &nval, &tval, &gval, &wval);
-               idx = NUM2INT(rb_Integer(ival));
-               if (idx < 0)
-                       idx = mol->npiatoms;
-               if (idx > mol->npiatoms || (func == 2 && idx == mol->npiatoms))
-                       rb_raise(rb_eMolbyError, "pi anchor index out of range");
-       }
-       np = StringValuePtr(nval);
-       strncpy(aname, np, 4);
-       tp = StringValuePtr(tval);
-       atype = AtomTypeEncodeToUInt(tp);
-       ig = IntGroupFromValue(gval);
-       if ((i = IntGroupGetCount(ig)) < 2)
-               rb_raise(rb_eMolbyError, "too few atoms are given (at least 2 are needed)");
-       ip = (Int *)malloc(sizeof(Int) * i);
-       for (i = 0; (j = IntGroupGetNthPoint(ig, i)) >= 0; i++) {
-               ip[i] = j;
-               if (j < 0 || j >= mol->natoms) {
-                       free(ip);
-                       rb_raise(rb_eMolbyError, "the atom index (%d) is out of range", j);
-               }
-       }
-       dp = (Double *)malloc(sizeof(Double) * i);
-       if (wval != Qnil) {
-               wval = rb_ary_to_ary(wval);
-               if (RARRAY_LEN(wval) < i) {
-                       free(ip);
-                       free(dp);
-                       rb_raise(rb_eMolbyError, "the number of floats is not sufficient (%d needed but only %d given)", i, RARRAY_LEN(wval));
-               }
-               for (j = 0; j < i; j++) {
-                       dp[j] = NUM2DBL(rb_Float(RARRAY_PTR(wval)[j]));
-               }
-       } else {
-               for (j = 0; j < i; j++) {
-                       dp[j] = 1.0 / i;
-               }
-       }
-       if (func == 0 || func == 1)  /*  Create or insert  */
-               MolActionCreateAndPerform(mol, gMolActionInsertOnePiAtom, idx, 4, aname, atype, i, ip, i, dp);
-       else  /*  Replace  */
-               MolActionCreateAndPerform(mol, gMolActionReplaceOnePiAtom, idx, 4, aname, atype, i, ip, i, dp);
-       free(ip);
-       free(dp);
-       return INT2NUM(idx);
-}
-
-static VALUE
-s_Molecule_CreatePiAnchor(int argc, VALUE *argv, VALUE self)
-{
-       return s_Molecule_CreateOrReplacePiAnchor(0, argc, argv, self);
-}
-
-static VALUE
-s_Molecule_InsertPiAnchor(int argc, VALUE *argv, VALUE self)
-{
-       return s_Molecule_CreateOrReplacePiAnchor(1, argc, argv, self);
-}
-
-static VALUE
-s_Molecule_ReplacePiAnchor(int argc, VALUE *argv, VALUE self)
-{
-       return s_Molecule_CreateOrReplacePiAnchor(2, argc, argv, self);
-}
-
-/*
- *  call-seq:
- *     remove_pi_anchor(idx)       -> Molecule
- *
- *  Remove the pi-anchor at the given index.
- */
-static VALUE
-s_Molecule_RemovePiAnchor(VALUE self, VALUE ival)
-{
-       Molecule *mol;
-       Int idx = NUM2INT(rb_Integer(ival));
-       Data_Get_Struct(self, Molecule, mol);
-       if (idx < 0 || idx >= mol->npiatoms)
-               rb_raise(rb_eMolbyError, "pi-anchor index (%d) is out of range (should be 0 <= index < %d)", idx, mol->npiatoms);
-       MolActionCreateAndPerform(mol, gMolActionRemoveOnePiAtom, idx);
-       return self;
-}
-
-/*
- *  call-seq:
- *     pi_anchor(idx) -> nil or [name, type, group, weights]
- *
- *  Return the attributes of the idx-th pi anchor if present, otherwise nil.
- */
-static VALUE
-s_Molecule_PiAnchorAtIndex(VALUE self, VALUE ival)
-{
-       Molecule *mol;
-       PiAtom *pp;
-       char aname[6], atype[8];
-       Int i, *connects;
-       VALUE cval, wval, rval;
-       Int idx = NUM2INT(rb_Integer(ival));
-    Data_Get_Struct(self, Molecule, mol);
-       if (idx < 0 || idx >= mol->npiatoms)
-               return Qnil;
-       pp = mol->piatoms + idx;
-       strncpy(aname, pp->aname, 4);
-       aname[4] = 0;
-       AtomTypeDecodeToString(pp->type, atype);
-       cval = rb_ary_new2(pp->connect.count);
-       connects = AtomConnectData(&pp->connect);
-       for (i = 0; i < pp->connect.count; i++)
-               rb_ary_store(cval, i, INT2NUM(connects[i]));
-       wval = rb_ary_new2(pp->ncoeffs);
-       for (i = 0; i < pp->ncoeffs; i++)
-               rb_ary_store(wval, i, rb_float_new(pp->coeffs[i]));
-       rval = rb_ary_new3(4, rb_str_new2(aname), rb_str_new2(atype), cval, wval);
-       return rval;
-}
-
-/*
- *  call-seq:
- *     pi_anchor_r(idx) -> Vector3D
- *
- *  Calculate the cartesian coordinate of the idx-th pi anchor if present.
- */
-static VALUE
-s_Molecule_PiAnchorR(VALUE self, VALUE ival)
-{
-       Molecule *mol;
-       Vector v;
-       Int idx = NUM2INT(rb_Integer(ival));
-    Data_Get_Struct(self, Molecule, mol);
-       if (idx < 0 || idx >= mol->npiatoms)
-               rb_raise(rb_eMolbyError, "no pi anchor present: %d", idx);
-       MoleculeCalculatePiAtomPosition(mol, idx);
-       v = mol->piatoms[idx].r;
-       return ValueFromVector(&v);
-}
-
-/*
- *  call-seq:
- *     pi_anchor_fract_r(idx) -> Vector3D
- *
- *  Calculate the fractional coordinate of the idx-th pi anchor if present.
- */
-static VALUE
-s_Molecule_PiAnchorFractR(VALUE self, VALUE ival)
-{
-       Molecule *mol;
-       Vector v;
-       Int idx = NUM2INT(rb_Integer(ival));
-    Data_Get_Struct(self, Molecule, mol);
-       if (idx < 0 || idx >= mol->npiatoms)
-               rb_raise(rb_eMolbyError, "no pi anchor present: %d", idx);
-       MoleculeCalculatePiAtomPosition(mol, idx);
-       v = mol->piatoms[idx].r;
-       if (mol->cell != NULL)
-               TransformVec(&v, mol->cell->rtr, &v);
-       return ValueFromVector(&v);
-}
-
-/*
- *  call-seq:
- *     count_pi_anchors -> Integer
- *
- *  Return the number of defined pi anchors.
- */
-static VALUE
-s_Molecule_CountPiAnchors(VALUE self)
-{
-       Molecule *mol;
-    Data_Get_Struct(self, Molecule, mol);
-       return INT2NUM(mol->npiatoms);
-}
-
-/*
- *  call-seq:
- *     create_pi_anchor_construct(n1, n2, n3 = nil, n4 = nil)       -> Integer
- *
- *  Create a bond, angle, or dihedral including one or more pi anchor points.
- *  The arguments can either be an atom representation (atom index, atom name, res_seq:name)
- *  or a pi-anchor representation (pi-anchor index, pi-anchor name).
- *  The pi-anchor index is represented by a negative integer -N, which corresponds to
- *  the (N-1)-th pi anchor.
- *  Returns the index for the newly created construct.
- *  This operation is undoable.
- */
-static VALUE
-s_Molecule_CreatePiAnchorConstruct(int argc, VALUE *argv, VALUE self)
-{
-    Molecule *mol;
-       VALUE vals[4];
-       Int ivals[4], i;
-       IntGroup *ig;
-    Data_Get_Struct(self, Molecule, mol);
-       rb_scan_args(argc, argv, "22", vals, vals + 1, vals + 2, vals + 3);
-       for (i = 0; i < 4; i++) {
-               if (vals[i] == Qnil)
-                       ivals[i] = -1;
-               else ivals[i] = s_Molecule_AtomOrPiAtomIndexFromValue(mol, vals[i]);
-       }
-       for (i = 0; i < 4; i++) {
-               if (ivals[i] >= ATOMS_MAX_NUMBER)
-                       break;
-       }
-       if (i == 4)
-               rb_raise(rb_eMolbyError, "No pi anchor is specified");
-       ig = IntGroupNewWithPoints(mol->npibonds, 1, -1);
-       MolActionCreateAndPerform(mol, gMolActionInsertPiBonds, ig, 4, ivals);
-       IntGroupRelease(ig);
-       return INT2NUM(mol->npibonds);
-}
-
-/*
- *  call-seq:
- *     remove_pi_anchor_constructs(IntGroup)       -> self
- *
- *  Remove pi anchor constructs (bond, angle, dihedral) with indices specified in IntGroup.
- *  This operation is undoable.
- */
-static VALUE
-s_Molecule_RemovePiAnchorConstructs(VALUE self, VALUE igval)
-{
-    Molecule *mol;
-       IntGroup *ig;
-    Data_Get_Struct(self, Molecule, mol);
-       ig = IntGroupFromValue(igval);
-       MolActionCreateAndPerform(mol, gMolActionRemovePiBonds, ig);
-       IntGroupRelease(ig);
-       return self;
-}
-
-/*
- *  call-seq:
- *     count_pi_anchor_constructs       -> Integer
- *
- *  Returns the number of pi anchor constructs.
- */
-static VALUE
-s_Molecule_CountPiAnchorConstructs(VALUE self)
-{
-    Molecule *mol;
-    Data_Get_Struct(self, Molecule, mol);
-       return INT2NUM(mol->npibonds);
-}
-
-/*
- *  call-seq:
- *     pi_anchor_construct(index)       -> Array of Integers
- *
- *  Returns the elements representing the index-th pi anchor constructs.
- */
-static VALUE
-s_Molecule_PiAnchorConstructAtIndex(VALUE self, VALUE ival)
-{
-    Molecule *mol;
-       Int idx, i, n;
-       VALUE vals[4];
-    Data_Get_Struct(self, Molecule, mol);
-       idx = NUM2INT(rb_Integer(ival));
-       if (idx < 0 || idx >= mol->npibonds)
-               rb_raise(rb_eMolbyError, "index of pi anchor constructs out of range (%d)", idx);
-       for (i = 0; i < 4; i++) {
-               n = mol->pibonds[idx * 4 + i];
-               if (n < 0)
-                       break;
-               if (n >= ATOMS_MAX_NUMBER)
-                       n = -(n - ATOMS_MAX_NUMBER) - 1;
-               vals[i] = INT2NUM(n);
-       }
-       return rb_ary_new4(i, vals);
-}
-#endif  /*  PIATOM  */
-
 /*
  *  call-seq:
  *     current       -> Molecule
@@ -10315,21 +9966,6 @@ Init_Molby(void)
        
        rb_define_method(rb_cMolecule, "create_pi_anchor", s_Molecule_CreatePiAnchor, -1);
 
-#if PIATOM
-       rb_define_method(rb_cMolecule, "replace_pi_anchor", s_Molecule_ReplacePiAnchor, -1);
-       rb_define_method(rb_cMolecule, "insert_pi_anchor", s_Molecule_InsertPiAnchor, -1);
-       rb_define_method(rb_cMolecule, "remove_pi_anchor", s_Molecule_RemovePiAnchor, 1);
-       rb_define_method(rb_cMolecule, "pi_anchor", s_Molecule_PiAnchorAtIndex, 1);
-       rb_define_method(rb_cMolecule, "pi_anchor_r", s_Molecule_PiAnchorR, 1);
-       rb_define_method(rb_cMolecule, "pi_anchor_fract_r", s_Molecule_PiAnchorFractR, 1);
-       rb_define_method(rb_cMolecule, "count_pi_anchors", s_Molecule_CountPiAnchors, 0);
-       rb_define_method(rb_cMolecule, "create_pi_anchor_construct", s_Molecule_CreatePiAnchorConstruct, -1);
-       rb_define_method(rb_cMolecule, "remove_pi_anchor_constructs", s_Molecule_RemovePiAnchorConstructs, 1);
-       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);