OSDN Git Service

Filter mode (development continued)
authortoshinagata1964 <toshinagata1964@a2be9bc6-48de-4e38-9406-05402d4bc13c>
Thu, 11 Oct 2012 15:45:00 +0000 (15:45 +0000)
committertoshinagata1964 <toshinagata1964@a2be9bc6-48de-4e38-9406-05402d4bc13c>
Thu, 11 Oct 2012 15:45:00 +0000 (15:45 +0000)
git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@291 a2be9bc6-48de-4e38-9406-05402d4bc13c

13 files changed:
MolLib/MD/MDCore.c
MolLib/MD/MDForce.c
MolLib/Molecule.c
MolLib/Molecule.h
MolLib/Ruby_bind/Molby_extern.h
MolLib/Ruby_bind/ruby_bind.c
wxSources/ConsoleFrame.cpp
wxSources/ConsoleFrame.h
wxSources/MoleculeView.cpp
wxSources/MyApp.cpp
wxSources/MyApp.h
wxSources/MyDocument.cpp
wxSources/MyDocument.h

index 336646b..8cf6340 100644 (file)
@@ -1203,7 +1203,7 @@ s_find_pibond_parameters(MDArena *arena)
                        if (ii >= ATOMS_MAX_NUMBER) {
                                ii -= ATOMS_MAX_NUMBER;
                                idx[j] = -1;
-                               if (ii < mol->npiatoms) {
+                               if (ii >= mol->npiatoms) {
                                        types[j] = -1;
                                        ptype = -1;
                                        break;
@@ -1743,6 +1743,7 @@ md_prepare(MDArena *arena, int check_only)
        s_find_angle_parameters(arena);
        s_find_dihedral_parameters(arena);
        s_find_improper_parameters(arena);
+       s_find_pibond_parameters(arena);
        
        if (arena->nmissing > 0) {
        /*      for (i = 0; i < nmissing; i++) {
@@ -3804,6 +3805,8 @@ md_arena_release(MDArena *arena)
                free(arena->vdw_par_i);
        if (arena->vdw_cache != NULL)
                free(arena->vdw_cache);
+       if (arena->pi_pars != NULL)
+               free(arena->pi_pars);
        if (arena->energies != NULL)
                free(arena->energies);
        if (arena->forces != NULL)
index f711b2a..89f363f 100644 (file)
@@ -96,6 +96,66 @@ s_calc_bond_force(MDArena *arena)
        }
 }
 
+static inline int
+s_calc_angle_force_one(const Vector *r1, const Vector *r2, const Vector *r3, const AnglePar *anp, Double *angle, Double *energy, Double *k1out, Vector *force1, Vector *force2, Vector *force3)
+{
+       Vector r21, r23, f1, f2, v1;
+       Double w1, w2, w3, cost, sint, t, k0, k1;
+       VecSub(r21, *r1, *r2);
+       VecSub(r23, *r3, *r2);
+       w1 = 1.0 / VecLength(r21);
+       w2 = 1.0 / VecLength(r23);
+       VecScaleSelf(r21, w1);
+       VecScaleSelf(r23, w2);
+       cost = VecDot(r21, r23);
+       if (cost > 0.999999 || cost < -0.999999) {
+               /*      printf("Cannot handle linear angle %d-%d-%d: skipped.\n", angles[0]+1, angles[1]+1, angles[2]+1); */
+               return -1;  /*  Cannot handle linear angle  */
+       }
+       sint = sqrt(1.0 - cost * cost);
+       /*      if (sint < 1e-5 && sint > -1e-5)
+        continue;  *//*  Cannot handle linear angle  */
+       t = atan2(sint, cost) - anp->a0;
+       k0 = anp->k * t * t;
+       
+       if (sint < 0.1 && sint > -0.1) {
+               /* ---- Method 1 ---- */
+               /* This is slower than method 2, but it is safer when sint is close to 0 */
+               k1 = -2 * anp->k * t;
+               VecCross(v1, r21, r23);
+               VecCross(f1, r21, v1);
+               w3 = w1 * k1 / VecLength(f1);
+               if (!isfinite(w3))
+                       return -1;
+               VecScaleSelf(f1, w3);
+               VecCross(f2, v1, r23);
+               w3 = w2 * k1 / VecLength(f2);
+               if (!isfinite(w3))
+                       return -1;
+               VecScaleSelf(f2, w3);
+       } else {
+               /* ---- Method 2 ---- */
+               k1 = -2 * anp->k * t / sint;
+               VecScale(f1, r21, cost);
+               VecDec(f1, r23);
+               w3 = w1 * k1;
+               VecScaleSelf(f1, w3);
+               VecScale(f2, r23, cost);
+               VecDec(f2, r21);
+               w3 = w2 * k1;
+               VecScaleSelf(f2, w3);
+       }
+       
+       *energy = k0;
+       *force1 = f1;
+       *force3 = f2;
+       *angle = t + anp->a0;
+       *k1out = k1;
+       VecInc(f1, f2);
+       VecScale(*force2, f1, -1);
+       return 0;
+}
+
 static void
 s_calc_angle_force(MDArena *arena)
 {
@@ -109,15 +169,18 @@ s_calc_angle_force(MDArena *arena)
        Vector *forces = &arena->forces[kAngleIndex * arena->mol->natoms];
        int i;
        for (i = 0; i < nangles; i++, angles += 3, angle_par_i++) {
-               Vector r21, r23, f1, f2, v1;
-               Double k0, k1, w1, w2, w3;
-               Double cost, sint, t;
+               Atom *ap1, *ap2, *ap3;
                AnglePar *anp;
+               Double en, ang, k1;
+               Vector f1, f2, f3;
                if (*angle_par_i < 0)
                        continue;  /*  Ignore this entry  */
+               ap1 = ATOM_AT_INDEX(ap, angles[0]);
+               ap2 = ATOM_AT_INDEX(ap, angles[1]);
+               ap3 = ATOM_AT_INDEX(ap, angles[2]);
        /*      if (angle_uniq != NULL && angle_uniq[i] >= 0)
                        continue; */ /*  Non-unique angle  */
-               if (ap[angles[0]].mm_exclude || ap[angles[1]].mm_exclude || ap[angles[2]].mm_exclude)
+               if (ap1->mm_exclude || ap2->mm_exclude || ap3->mm_exclude)
                        continue;  /*  Skip non-occupied atoms  */
                if (arena->nalchem_flags > 0) {
                        if (angles[0] < arena->nalchem_flags && angles[1] < arena->nalchem_flags
@@ -126,62 +189,124 @@ s_calc_angle_force(MDArena *arena)
                                continue;  /*  Interaction between vanishing and appearing groups is ignored  */
                }
                anp = angle_pars + *angle_par_i;
-               VecSub(r21, ap[angles[0]].r, ap[angles[1]].r);
-               VecSub(r23, ap[angles[2]].r, ap[angles[1]].r);
-               w1 = 1.0 / VecLength(r21);
-               w2 = 1.0 / VecLength(r23);
-               VecScaleSelf(r21, w1);
-               VecScaleSelf(r23, w2);
-               cost = VecDot(r21, r23);
-               if (cost > 0.999999 || cost < -0.999999) {
-               /*      printf("Cannot handle linear angle %d-%d-%d: skipped.\n", angles[0]+1, angles[1]+1, angles[2]+1); */
-                       continue;  /*  Cannot handle linear angle  */
-               }
-               sint = sqrt(1.0 - cost * cost);
-       /*      if (sint < 1e-5 && sint > -1e-5)
-                       continue;  *//*  Cannot handle linear angle  */
-               t = atan2(sint, cost) - anp->a0;
-               k0 = anp->k * t * t;
-       
-               if (sint < 0.1 && sint > -0.1) {
-                       /* ---- Method 1 ---- */
-                       /* This is slower than method 2, but it is safer when sint is close to 0 */
-                       k1 = -2 * anp->k * t;
-                       VecCross(v1, r21, r23);
-                       VecCross(f1, r21, v1);
-                       w3 = w1 * k1 / VecLength(f1);
-                       if (!isfinite(w3))
-                               continue;
-                       VecScaleSelf(f1, w3);
-                       VecCross(f2, v1, r23);
-                       w3 = w2 * k1 / VecLength(f2);
-                       if (!isfinite(w3))
-                               continue;
-                       VecScaleSelf(f2, w3);
-               } else {
-                       /* ---- Method 2 ---- */
-                       k1 = -2 * anp->k * t / sint;
-                       VecScale(f1, r21, cost);
-                       VecDec(f1, r23);
-                       w3 = w1 * k1;
-                       VecScaleSelf(f1, w3);
-                       VecScale(f2, r23, cost);
-                       VecDec(f2, r21);
-                       w3 = w2 * k1;
-                       VecScaleSelf(f2, w3);
-               }
-
-               *energies += k0;
+               if (s_calc_angle_force_one(&(ap1->r), &(ap2->r), &(ap3->r), angle_pars + *angle_par_i, &ang, &en, &k1, &f1, &f2, &f3) != 0)
+                       continue;
+               *energies += en;
                VecInc(forces[angles[0]], f1);
-               VecInc(forces[angles[2]], f2);
-               VecInc(f1, f2);
-               VecDec(forces[angles[1]], f1);
-
+               VecInc(forces[angles[1]], f2);
+               VecInc(forces[angles[2]], f3);
                if (arena->debug_result && arena->debug_output_level > 1) {
-                       fprintf(arena->debug_result, "angle force %d-%d-%d: a=%f, a0=%f, k1=%f, {%f %f %f}, {%f %f %f}\n", angles[0]+1, angles[1]+1, angles[2]+1, (t+anp->a0)*180/PI, anp->a0*180/PI, k1, f1.x, f1.y, f1.z, f2.x, f2.y, f2.z);
+                       fprintf(arena->debug_result, "angle force %d-%d-%d: a=%f, a0=%f, k1=%f, {%f %f %f}, {%f %f %f}\n", angles[0]+1, angles[1]+1, angles[2]+1, ang*180/PI, anp->a0*180/PI, k1, -f2.x, -f2.y, -f2.z, f3.x, f3.y, f3.z);
                }
+       }
+}
 
+static inline int
+s_calc_dihedral_force_one(const Vector *r1, const Vector *r2, const Vector *r3, const Vector *r4, const TorsionPar *tp, Double *phi, Double *energy, Double *k1out, Vector *force1, Vector *force2, Vector *force3, Vector *force4)
+{
+       Vector r21, r32, r43;
+       Vector v1, v2, v3;
+       Double w1, w2, w3, k0, k1;
+       Double cosphi, sinphi;
+       Vector f1, f2, f3;
+       Int n;
+       VecSub(r21, *r1, *r2);
+       VecSub(r32, *r2, *r3);
+       VecSub(r43, *r3, *r4);
+       VecCross(v1, r21, r32);
+       VecCross(v2, r32, r43);
+       VecCross(v3, r32, v1);
+       w1 = VecLength(v1);
+       w2 = VecLength(v2);
+       w3 = VecLength(v3);
+       if (w1 < 1e-5 || w2 < 1e-5 || w3 < 1e-5)
+               return -1;  /*  The dihedral cannot be defined  */
+       w1 = 1.0/w1;
+       w2 = 1.0/w2;
+       w3 = 1.0/w3;
+       VecScaleSelf(v1, w1);
+       VecScaleSelf(v2, w2);
+       VecScaleSelf(v3, w3);
+       
+       /*  The dihedral angle value  */
+       cosphi = VecDot(v1, v2);
+       sinphi = VecDot(v3, v2);
+       *phi = -atan2(sinphi, cosphi);
+       
+       /*  Repeat for multiple dihedral terms  */
+       k0 = k1 = 0.0;
+       for (n = 0; n < tp->mult; n++) {
+               Double k = tp->k[n];
+               Double phi0 = tp->phi0[n];
+               Int period = tp->period[n];
+               if (period > 0) {
+                       k0 += k * (1 + cos(period * *phi + phi0));
+                       k1 -= period * k * sin(period * *phi + phi0);
+               } else {
+                       /*  A simple quadratic term  */
+                       phi0 = *phi - phi0;
+                       if (phi0 < -PI)
+                               phi0 += 2 * PI;
+                       else if (phi0 > PI)
+                               phi0 -= 2 * PI;
+                       k0 += k * phi0 * phi0;
+                       k1 += 2 * k * phi0;
+               }
        }
+       
+       if (sinphi < -0.1 || sinphi > 0.1) {
+               /*  The sin form  */
+               Vector v4, v5, vw0;
+               VecScale(v4, v1, cosphi);
+               VecDec(v4, v2);
+               VecScaleSelf(v4, w1);
+               VecScale(v5, v2, cosphi);
+               VecDec(v5, v1);
+               VecScaleSelf(v5, w2);
+               k1 /= sinphi;
+               VecCross(f1, r32, v4);
+               VecCross(f3, v5, r32);
+               VecCross(f2, v4, r21);
+               VecCross(vw0, r43, v5);
+               VecInc(f2, vw0);
+               VecScaleSelf(f1, k1);
+               VecScaleSelf(f2, k1);
+               VecScaleSelf(f3, k1);
+       } else {
+               /*  The cos form  */
+               Vector v6, v7, vw1, vw2;
+               VecScale(v6, v3, sinphi);
+               VecDec(v6, v2);
+               VecScaleSelf(v6, w3);
+               VecScale(v7, v2, sinphi);
+               VecDec(v7, v3);
+               VecScaleSelf(v7, w2);
+               k1 = -k1 / cosphi;
+               VecCross(vw1, v6, r32);
+               VecCross(f1, r32, vw1);
+               VecCross(f3, v7, r32);
+               VecCross(f2, r43, v7);
+               VecCross(vw1, r32, r21);
+               VecCross(vw2, v6, vw1);
+               VecInc(f2, vw2);
+               VecCross(vw1, r32, v6);
+               VecCross(vw2, r21, vw1);
+               VecInc(f2, vw2);
+               VecScaleSelf(f1, k1);
+               VecScaleSelf(f2, k1);
+               VecScaleSelf(f3, k1);
+       }
+       /*      if (dihedral_uniq != NULL)
+        k0 *= -dihedral_uniq[i]; */
+       *energy = k0;
+       *k1out = k1;
+       *force1 = f1;
+       VecDec(f1, f2);
+       VecScale(*force2, f1, -1);
+       VecDec(f2, f3);
+       VecScale(*force3, f2, -1);
+       VecScale(*force4, f3, -1);
+       return 0;
 }
 
 static void
@@ -191,19 +316,20 @@ s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedra
        if (arena->mol->nsyms == 0)
                dihedral_uniq = NULL;  /*  Ignore the symmetry info  */
        for (i = 0; i < ndihedrals; i++, dihedrals += 4, dihedral_par_i++) {
-               Vector r21, r32, r43;
-               Vector v1, v2, v3;
-               Double w1, w2, w3, k0, k1;
-               Double cosphi, sinphi, phi;
-               Vector f1, f2, f3;
+               Atom *ap1, *ap2, *ap3, *ap4;
+               Double phi, en, k1;
+               Vector f1, f2, f3, f4;
                TorsionPar *tp;
-               int n;
 
                if (*dihedral_par_i < 0)
                        continue;  /*  Ignore this entry  */
        /*      if (dihedral_uniq != NULL && dihedral_uniq[i] >= 0)
                        continue;  *//*  Non-unique dihedral  */
-               if (ap[dihedrals[0]].mm_exclude || ap[dihedrals[1]].mm_exclude || ap[dihedrals[2]].mm_exclude || ap[dihedrals[3]].mm_exclude)
+               ap1 = ATOM_AT_INDEX(ap, dihedrals[0]);
+               ap2 = ATOM_AT_INDEX(ap, dihedrals[1]);
+               ap3 = ATOM_AT_INDEX(ap, dihedrals[2]);
+               ap4 = ATOM_AT_INDEX(ap, dihedrals[3]);
+               if (ap1->mm_exclude || ap2->mm_exclude || ap3->mm_exclude || ap4->mm_exclude)
                        continue;  /*  Skip non-occupied atoms  */
                if (arena->nalchem_flags > 0) {
                        if (dihedrals[0] < arena->nalchem_flags && dihedrals[1] < arena->nalchem_flags
@@ -213,6 +339,10 @@ s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedra
                                continue;  /*  Interaction between vanishing and appearing groups is ignored  */
                }
                tp = dihedral_pars + *dihedral_par_i;
+               if (s_calc_dihedral_force_one(&(ap1->r), &(ap2->r), &(ap3->r), &(ap4->r), tp, &phi, &en, &k1, &f1, &f2, &f3, &f4) != 0)
+                       continue;
+               
+               /*
                VecSub(r21, ap[dihedrals[0]].r, ap[dihedrals[1]].r);
                VecSub(r32, ap[dihedrals[1]].r, ap[dihedrals[2]].r);
                VecSub(r43, ap[dihedrals[2]].r, ap[dihedrals[3]].r);
@@ -223,7 +353,7 @@ s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedra
                w2 = VecLength(v2);
                w3 = VecLength(v3);
                if (w1 < 1e-5 || w2 < 1e-5 || w3 < 1e-5)
-                       continue;  /*  The dihedral cannot be defined  */
+                       continue;  //  The dihedral cannot be defined 
                w1 = 1.0/w1;
                w2 = 1.0/w2;
                w3 = 1.0/w3;
@@ -231,12 +361,12 @@ s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedra
                VecScaleSelf(v2, w2);
                VecScaleSelf(v3, w3);
                
-               /*  The dihedral angle value  */
+               //  The dihedral angle value
                cosphi = VecDot(v1, v2);
                sinphi = VecDot(v3, v2);
                phi = -atan2(sinphi, cosphi);
                
-               /*  Repeat for multiple dihedral terms  */
+               //  Repeat for multiple dihedral terms
                k0 = k1 = 0.0;
                for (n = 0; n < tp->mult; n++) {
                        Double k = tp->k[n];
@@ -246,7 +376,7 @@ s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedra
                                k0 += k * (1 + cos(period * phi + phi0));
                                k1 -= period * k * sin(period * phi + phi0);
                        } else {
-                               /*  A simple quadratic term  */
+                               //  A simple quadratic term
                                phi -= phi0;
                                if (phi < -PI)
                                        phi += 2 * PI;
@@ -258,7 +388,7 @@ s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedra
                }
                
                if (sinphi < -0.1 || sinphi > 0.1) {
-                       /*  The sin form  */
+                       //  The sin form
                        Vector v4, v5, vw0;
                        VecScale(v4, v1, cosphi);
                        VecDec(v4, v2);
@@ -276,7 +406,7 @@ s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedra
                        VecScaleSelf(f2, k1);
                        VecScaleSelf(f3, k1);
                } else {
-                       /*  The cos form  */
+                       //  The cos form
                        Vector v6, v7, vw1, vw2;
                        VecScale(v6, v3, sinphi);
                        VecDec(v6, v2);
@@ -299,8 +429,6 @@ s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedra
                        VecScaleSelf(f2, k1);
                        VecScaleSelf(f3, k1);
                }
-       /*      if (dihedral_uniq != NULL)
-                       k0 *= -dihedral_uniq[i]; */
                *energies += k0;
                VecInc(forces[dihedrals[0]], f1);
                VecDec(f1, f2);
@@ -308,7 +436,12 @@ s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedra
                VecDec(f2, f3);
                VecDec(forces[dihedrals[2]], f2);
                VecDec(forces[dihedrals[3]], f3);
-
+*/
+               *energies += en;
+               VecInc(forces[dihedrals[0]], f1);
+               VecInc(forces[dihedrals[1]], f2);
+               VecInc(forces[dihedrals[2]], f3);
+               VecInc(forces[dihedrals[3]], f4);
                if (arena->debug_result && arena->debug_output_level > 1) {
                        fprintf(arena->debug_result, "dihedral(improper) force %d-%d-%d-%d: phi=%f, k1=%f, {%f %f %f}, {%f %f %f}, {%f %f %f}\n", dihedrals[0]+1, dihedrals[1]+1, dihedrals[2]+1, dihedrals[3]+1, phi*180/PI, k1, f1.x, f1.y, f1.z, f2.x, f2.y, f2.z, f3.x, f3.y, f3.z);
                }
@@ -335,6 +468,67 @@ 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);
 }
 
+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;
+//     Double *energies = &arena->energies[kAngleIndex];
+//     Vector *forces = &arena->forces[kAngleIndex * arena->mol->natoms];
+       Int i;
+       for (i = 0; i < npibonds; i++, pibonds += 4, pars++) {
+               Vector r[4];
+               Vector f[4];
+               PiAtom *pp[4];
+               Int j, 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;
+                               MoleculeCalculatePiAtomPosition(arena->mol, n, &r[j]);
+                               pp[j] = arena->mol->piatoms + n;
+                       }
+               }
+               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);
+                       }
+               } 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);
+                       }
+               } 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);
+                       }
+               } else continue;
+       }
+       
+}
+
 /*  ==================================================================== */
 /*  md_check_verlet_list: only for debugging the verlet list generation  */
 /*  ==================================================================== */
index 7921718..531301a 100755 (executable)
@@ -10050,6 +10050,25 @@ MoleculeFlushFrames(Molecule *mp)
        return nframes;
 }
 
+#pragma mark ====== Pi Atoms ======
+
+int
+MoleculeCalculatePiAtomPosition(Molecule *mol, int idx, Vector *vp)
+{
+       PiAtom *pp;
+       Int i, *cp;
+       if (mol == NULL || idx < 0 || idx >= mol->npiatoms)
+               return -1;
+       pp = mol->piatoms + idx;
+       cp = AtomConnectData(&pp->connect);
+       VecZero(*vp);
+       for (i = pp->connect.count - 1; i >= 0; i--) {
+               Vector r = ATOM_AT_INDEX(mol->atoms, cp[i])->r;
+               VecScaleInc(*vp, r, pp->coeffs[i]);
+       }
+       return idx;
+}
+
 #pragma mark ====== MO calculation ======
 
 /*  Calculate an MO value for a single point.  */
index 9f14e4d..678465e 100755 (executable)
@@ -513,6 +513,8 @@ int MoleculeRemoveFrames(Molecule *mp, IntGroup *group, Vector *outFrame, Vector
 int MoleculeSelectFrame(Molecule *mp, int frame, int copyback);
 int MoleculeFlushFrames(Molecule *mp);
 
+int MoleculeCalculatePiAtomPosition(Molecule *mol, int idx, Vector *vp);
+
 int MoleculeCalcMO(Molecule *mp, Int mono, const Vector *op, const Vector *dxp, const Vector *dyp, const Vector *dzp, Int nx, Int ny, Int nz, int (*callback)(double progress, void *ref), void *ref);
 int MoleculeGetDefaultMOGrid(Molecule *mp, Int npoints, Vector *op, Vector *xp, Vector *yp, Vector *zp, Int *nx, Int *ny, Int *nz);
 const Cube *MoleculeGetCubeAtIndex(Molecule *mp, Int index);
index e226200..4b6e9ab 100755 (executable)
@@ -62,7 +62,8 @@ STUB RubyValue MyAppCallback_executeScriptFromFile(const char *path, int *status
 STUB int MyAppCallback_callSubProcess(const char *cmdline, const char *procname);
 STUB void MyAppCallback_beginUndoGrouping(void);
 STUB void MyAppCallback_endUndoGrouping(void);
-
+STUB int MyAppCallback_switchToFilterMode(void);
+       
 #ifdef __cplusplus
 }
 #endif
index 465604a..d9ee34c 100644 (file)
@@ -684,6 +684,31 @@ s_Kernel_CallSubProcess(VALUE self, VALUE cmd, VALUE procname)
        return INT2NUM(n);
 }
 
+/*
+ *  call-seq:
+ *     enable_filter_mode -> Integer
+ *
+ *  Switch to a filter mode. Should be invoked from within a script file.
+ *  Returns true when successfully switched, false when already in filter mode.
+ *  No Molecule should be open as a view; otherwise, an exception is raised.
+ */
+static VALUE
+s_Kernel_EnableFilterMode(VALUE self)
+{
+       int n = MyAppCallback_switchToFilterMode();
+       if (n == 0)
+               return Qtrue;
+       else if (n == 1)
+               return Qfalse;
+       else if (n == -1)
+               rb_raise(rb_eMolbyError, "To switch to filter mode, all molecule should be closed");
+       else if (n == -2)
+               rb_raise(rb_eMolbyError, "No script file is specified for filter mode");
+       else
+               rb_raise(rb_eMolbyError, "Cannot switch to filter mode");
+       return Qnil;  /*  Dummy to keep compiler happy  */
+}
+
 #pragma mark ====== User defaults ======
 
 /*
@@ -4241,6 +4266,7 @@ s_Molecule_AtomOrPiAtomIndexFromValue(Molecule *mol, VALUE val)
                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));
        }
@@ -4265,8 +4291,8 @@ s_Molecule_AtomIndexFromValue(Molecule *mol, VALUE val)
 {
        int n = s_Molecule_AtomOrPiAtomIndexFromValue(mol, val);
        if (n >= ATOMS_MAX_NUMBER) {
-               char *p = StringValuePtr(val);
-               rb_raise(rb_eMolbyError, "no such atom: %s", p);
+               val = rb_inspect(val);
+               rb_raise(rb_eMolbyError, "no such atom: %s", StringValuePtr(val));
        }
        return n;
 }
@@ -9508,6 +9534,46 @@ s_Molecule_PiAnchorAtIndex(VALUE self, VALUE ival)
 
 /*
  *  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);
+       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);
+       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.
@@ -9954,6 +10020,8 @@ Init_Molby(void)
        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);
@@ -10079,6 +10147,7 @@ Init_Molby(void)
        rb_define_method(rb_mKernel, "call_subprocess", s_Kernel_CallSubProcess, 2);
        rb_define_method(rb_mKernel, "message_box", s_Kernel_MessageBox, -1);
        rb_define_method(rb_mKernel, "error_message_box", s_Kernel_ErrorMessageBox, 1);
+       rb_define_method(rb_mKernel, "enable_filter_mode", s_Kernel_EnableFilterMode, 0);
        
        s_ID_equal = rb_intern("==");
 }
index 5e2675f..7f49aa6 100755 (executable)
@@ -62,7 +62,7 @@ ConsoleFrame::OnCreate()
        textCtrl = new wxTextCtrl(this, wxID_ANY, _T(""), wxPoint(0, 0), wxSize(100, 100), wxTE_MULTILINE | wxTE_RICH);
 
        wxBoxSizer *consoleSizer = new wxBoxSizer(wxHORIZONTAL);
-       consoleSizer->Add(textCtrl, 1, wxEXPAND | wxALL, 2);
+       consoleSizer->Add(textCtrl, 1, wxEXPAND);
        this->SetSizer(consoleSizer);
        consoleSizer->SetSizeHints(this);
        
@@ -308,8 +308,9 @@ ConsoleFrame::OnKeyDown(wxKeyEvent &event)
 }
 
 void
-ConsoleFrame::EmptyBuffer()
+ConsoleFrame::EmptyBuffer(bool showRubyPrompt)
 {
        textCtrl->Clear();
-       MyAppCallback_showRubyPrompt();
+       if (showRubyPrompt)
+               MyAppCallback_showRubyPrompt();
 }
index da9c3df..e78fd2a 100755 (executable)
@@ -48,7 +48,7 @@ public:
        void OnCopy(wxCommandEvent &event);
        void OnPaste(wxCommandEvent &event);
        void OnClear(wxCommandEvent &event);
-       void EmptyBuffer();
+       void EmptyBuffer(bool showRubyPrompt = true);
 
 private:
        DECLARE_EVENT_TABLE()
index 5192a54..707eb7e 100755 (executable)
@@ -95,6 +95,7 @@ MoleculeView::OnCreate(wxDocument *doc, long WXUNUSED(flags) )
 {
        int i;
 
+       
        // Make a document frame
        frame = new wxDocMDIChildFrame(doc, this, GetMainFrame(), wxID_ANY, _T("New Molby Document"),
                                                   wxPoint(10, 24), wxSize(680, 400),
index f5da639..91e01a3 100755 (executable)
@@ -74,6 +74,8 @@
 
 static char *sLastBuildString = "";
 
+static const char *sExecutingRubyScriptPath = NULL;
+
 MyFrame *frame = (MyFrame *) NULL;
 
 IMPLEMENT_APP(MyApp)
@@ -658,7 +660,8 @@ MyApp::MacHandleAEODoc(const WXEVENTREF event, WXEVENTREF WXUNUSED(reply))
 
 #endif
 
-int MyApp::OnExit(void)
+int
+MyApp::OnExit(void)
 {
        SaveDefaultSettings();
     delete m_docManager;
@@ -669,6 +672,111 @@ int MyApp::OnExit(void)
     return 0;
 }
 
+static void
+sModifyMenuForFilterMode(wxMenuBar *mbar)
+{
+       int idx, i, n, id;
+       wxMenu *menu;
+       wxMenuItem *item;
+       idx = mbar->FindMenu(wxT("Show"));
+       if (idx != wxNOT_FOUND)
+               delete mbar->Remove(idx);
+       idx = mbar->FindMenu(wxT("MM/MD"));
+       if (idx != wxNOT_FOUND)
+               delete mbar->Remove(idx);
+       idx = mbar->FindMenu(wxT("QChem"));
+       if (idx != wxNOT_FOUND)
+               delete mbar->Remove(idx);
+       idx = mbar->FindMenu(wxT("Script"));
+       if (idx != wxNOT_FOUND) {
+               menu = mbar->GetMenu(idx);
+               n = menu->GetMenuItemCount();
+               for (i = n - 1; i >= 0; i--) {
+                       item = menu->FindItemByPosition(i);
+                       id = item->GetId();
+                       if (id != myMenuID_OpenConsoleWindow && id != myMenuID_EmptyConsoleWindow && id != myMenuID_ExecuteScript) {
+                               menu->Remove(item);
+                               delete item;
+                       }
+               }
+       }
+       
+       idx = mbar->FindMenu(wxT("Edit"));
+       if (idx != wxNOT_FOUND) {
+               menu = mbar->GetMenu(idx);
+               n = menu->GetMenuItemCount();
+               for (i = n - 1; i >= 0; i--) {
+                       item = menu->FindItemByPosition(i);
+                       id = item->GetId();
+                       if (id == wxID_SELECTALL)
+                               break;
+                       menu->Remove(item);
+                       delete item;
+               }
+       }
+       
+       idx = mbar->FindMenu(wxT("File"));
+       if (idx != wxNOT_FOUND) {
+               menu = mbar->GetMenu(idx);
+               n = menu->GetMenuItemCount();
+               for (i = n - 1; i >= 0; i--) {
+                       item = menu->FindItemByPosition(i);
+                       id = item->GetId();
+                       if (id != wxID_OPEN && id != wxID_EXIT) {
+                               menu->Remove(item);
+                               delete item;
+                       }
+               }
+       }
+       
+}
+
+int
+MyApp::SwitchToFilterMode(const char *filterScriptName)
+{
+       if (IsFilterMode())
+               return 1;   /*  Already filter mode  */
+       if (m_docManager->GetCurrentView() != NULL)
+               return -1;  /*  Molecule is open: cannot switch  */
+
+       /*  Remove menu items except for absolutely necessary  */
+       sModifyMenuForFilterMode(GetMainFrame()->GetMenuBar());
+       sModifyMenuForFilterMode(consoleFrame->GetMenuBar());
+       
+       /*  Record the name of the filter script  */
+       char *p;
+       m_filterScriptName = strdup(filterScriptName);
+       p = strrchr(m_filterScriptName, '/');
+#if defined(__WXMSW__)
+       if (p == NULL)
+               p = strrchr(m_filterScriptName, '\\');
+#endif
+       if (p == NULL)
+               p = m_filterScriptName;
+       else p++;
+       m_filterScriptBaseName = strdup(p);
+
+       /*  Resize the console window and show startup message  */
+       int width, height;
+       consoleFrame->GetClientSize(&width, &height);
+       if (width < 640)
+               width = 640;
+       if (height < 480)
+               height = 480;
+       consoleFrame->EmptyBuffer(false);
+       consoleFrame->SetClientSize(width, height);
+       MyAppCallback_setConsoleColor(3);
+       MyAppCallback_showScriptMessage("***********************************************************\n");
+       MyAppCallback_showScriptMessage(" Molby is now running in the filter mode.\n");
+       MyAppCallback_showScriptMessage(" When you open files, they will be processed by\n");
+       MyAppCallback_showScriptMessage(" the script '%s'.\n", m_filterScriptBaseName);
+       MyAppCallback_showScriptMessage(" You can process any file, not necessarily molecular files.\n");
+       MyAppCallback_showScriptMessage("***********************************************************\n");
+       MyAppCallback_setConsoleColor(0);
+       
+       return 0;
+}
+
 int
 MyApp::AppendConsoleMessage(const char *mes)
 {
@@ -1159,7 +1267,7 @@ MyApp::OnOpenFilesByIPC(wxCommandEvent& event)
 }
 
 bool
-MyApp::OnOpenFiles(wxString &files)
+MyApp::OnOpenFiles(const wxString &files)
 {
        Int start, end;
        Int nargs = 0;
@@ -1188,9 +1296,10 @@ MyApp::OnOpenFiles(wxString &files)
                                                Molby_showError(status);
                                        return false;
                                }
+                       } else {
+                               if (NULL == wxGetApp().DocManager()->CreateDocument(file, wxDOC_SILENT))
+                                       success = false;
                        }
-                       if (NULL == wxGetApp().DocManager()->CreateDocument(file, wxDOC_SILENT))
-                               success = false;
                }
                if (end == wxString::npos)
                        break;
@@ -1769,6 +1878,7 @@ RubyValue
 MyAppCallback_executeScriptFromFile(const char *cpath, int *status)
 {
        RubyValue retval;
+       const char *old_cpath;
        wxString cwd = wxFileName::GetCwd();
        wxString path(cpath, wxConvFile);
        char *p = strdup(cpath);
@@ -1844,7 +1954,11 @@ MyAppCallback_executeScriptFromFile(const char *cpath, int *status)
                }
        }
        
+       old_cpath = sExecutingRubyScriptPath;
+       sExecutingRubyScriptPath = cpath;
        retval = Molby_evalRubyScriptOnMolecule(script, MoleculeCallback_currentMolecule(), pp, status);
+       sExecutingRubyScriptPath = old_cpath;
+       
        free(script);
        free(p);
        wxFileName::SetCwd(cwd);
@@ -1873,3 +1987,10 @@ int MyAppCallback_callSubProcess(const char *cmdline, const char *procname)
 {
        return wxGetApp().CallSubProcess(cmdline, procname);
 }
+
+int MyAppCallback_switchToFilterMode(void)
+{
+       if (sExecutingRubyScriptPath == NULL)
+               return -2;  /*  Not invoked from file  */
+       return wxGetApp().SwitchToFilterMode(sExecutingRubyScriptPath);
+}
index d0ba1ef..fec47f3 100755 (executable)
@@ -180,8 +180,9 @@ class MyApp: public wxApp
        void RequestOpenFilesByIPC(wxString& files);
        void OnOpenFilesByIPC(wxCommandEvent& event);
        
-       bool OnOpenFiles(wxString &files);
+       bool OnOpenFiles(const wxString &files);
 
+       int SwitchToFilterMode(const char *filterScriptName);
        bool IsFilterMode() { return (m_filterScriptName != NULL); }
        const char *FilterScriptBaseName() { return m_filterScriptBaseName; }
 
@@ -218,9 +219,9 @@ protected:
 
        wxString *m_pendingFilesToOpen;  /*  Files to be processed by OnOpenFilesByIPC()  */
 
-       char *m_filterScriptName;         /*  Ruby script when invoked as a filter mode  */
-       char *m_filterScriptBaseName;     /*  The file name (without directories) of the filter script  */
-
+       char *m_filterScriptName;        /*  Ruby script when invoked as a filter mode  */
+       char *m_filterScriptBaseName;    /*  The file name (without directories) of the filter script  */
+       
 #if defined(__WXMSW__)
 public:
        wxSingleInstanceChecker *m_checker;
index 896dc00..7bc7744 100755 (executable)
@@ -262,6 +262,18 @@ MyDocument::DoOpenDocument(const wxString& file)
        return true;
 }
 
+/*  Override to intercept view creation for running script  */
+bool
+MyDocument::OnCreate(const wxString& path, long flags)
+{
+       if (path.EndsWith(wxT(".rb")) || path.EndsWith(wxT(".mrb"))) {
+               wxGetApp().OnOpenFiles(path);
+               return false;  /*  This document will be deleted  */
+       } else {
+               return wxDocument::OnCreate(path, flags);
+       }
+}
+
 void
 MyDocument::OnImport(wxCommandEvent& event)
 {
index 4d06cef..af52278 100755 (executable)
@@ -133,6 +133,7 @@ public:
  protected:
        virtual bool DoSaveDocument(const wxString& file);
        virtual bool DoOpenDocument(const wxString& file);
+       virtual bool OnCreate(const wxString& path, long flags);
 
  private:
        DECLARE_DYNAMIC_CLASS(MyDocument)