From dd4317536b287f624d1e894276aa28da185602b7 Mon Sep 17 00:00:00 2001 From: toshinagata1964 Date: Thu, 11 Oct 2012 15:45:00 +0000 Subject: [PATCH] Filter mode (development continued) git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@291 a2be9bc6-48de-4e38-9406-05402d4bc13c --- MolLib/MD/MDCore.c | 5 +- MolLib/MD/MDForce.c | 336 +++++++++++++++++++++++++++++++--------- MolLib/Molecule.c | 19 +++ MolLib/Molecule.h | 2 + MolLib/Ruby_bind/Molby_extern.h | 3 +- MolLib/Ruby_bind/ruby_bind.c | 73 ++++++++- wxSources/ConsoleFrame.cpp | 7 +- wxSources/ConsoleFrame.h | 2 +- wxSources/MoleculeView.cpp | 1 + wxSources/MyApp.cpp | 129 ++++++++++++++- wxSources/MyApp.h | 9 +- wxSources/MyDocument.cpp | 12 ++ wxSources/MyDocument.h | 1 + 13 files changed, 512 insertions(+), 87 deletions(-) diff --git a/MolLib/MD/MDCore.c b/MolLib/MD/MDCore.c index 336646b..8cf6340 100644 --- a/MolLib/MD/MDCore.c +++ b/MolLib/MD/MDCore.c @@ -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) diff --git a/MolLib/MD/MDForce.c b/MolLib/MD/MDForce.c index f711b2a..89f363f 100644 --- a/MolLib/MD/MDForce.c +++ b/MolLib/MD/MDForce.c @@ -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 */ /* ==================================================================== */ diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c index 7921718..531301a 100755 --- a/MolLib/Molecule.c +++ b/MolLib/Molecule.c @@ -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. */ diff --git a/MolLib/Molecule.h b/MolLib/Molecule.h index 9f14e4d..678465e 100755 --- a/MolLib/Molecule.h +++ b/MolLib/Molecule.h @@ -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); diff --git a/MolLib/Ruby_bind/Molby_extern.h b/MolLib/Ruby_bind/Molby_extern.h index e226200..4b6e9ab 100755 --- a/MolLib/Ruby_bind/Molby_extern.h +++ b/MolLib/Ruby_bind/Molby_extern.h @@ -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 diff --git a/MolLib/Ruby_bind/ruby_bind.c b/MolLib/Ruby_bind/ruby_bind.c index 465604a..d9ee34c 100644 --- a/MolLib/Ruby_bind/ruby_bind.c +++ b/MolLib/Ruby_bind/ruby_bind.c @@ -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("=="); } diff --git a/wxSources/ConsoleFrame.cpp b/wxSources/ConsoleFrame.cpp index 5e2675f..7f49aa6 100755 --- a/wxSources/ConsoleFrame.cpp +++ b/wxSources/ConsoleFrame.cpp @@ -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(); } diff --git a/wxSources/ConsoleFrame.h b/wxSources/ConsoleFrame.h index da9c3df..e78fd2a 100755 --- a/wxSources/ConsoleFrame.h +++ b/wxSources/ConsoleFrame.h @@ -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() diff --git a/wxSources/MoleculeView.cpp b/wxSources/MoleculeView.cpp index 5192a54..707eb7e 100755 --- a/wxSources/MoleculeView.cpp +++ b/wxSources/MoleculeView.cpp @@ -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), diff --git a/wxSources/MyApp.cpp b/wxSources/MyApp.cpp index f5da639..91e01a3 100755 --- a/wxSources/MyApp.cpp +++ b/wxSources/MyApp.cpp @@ -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); +} diff --git a/wxSources/MyApp.h b/wxSources/MyApp.h index d0ba1ef..fec47f3 100755 --- a/wxSources/MyApp.h +++ b/wxSources/MyApp.h @@ -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; diff --git a/wxSources/MyDocument.cpp b/wxSources/MyDocument.cpp index 896dc00..7bc7744 100755 --- a/wxSources/MyDocument.cpp +++ b/wxSources/MyDocument.cpp @@ -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) { diff --git a/wxSources/MyDocument.h b/wxSources/MyDocument.h index 4d06cef..af52278 100755 --- a/wxSources/MyDocument.h +++ b/wxSources/MyDocument.h @@ -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) -- 2.11.0