md_center_of_mass(arena, &(arena->initial_center));
}
+/* Set the alchemical flags */
+/* Independent with arena->xmol, mol */
+int
+md_set_alchemical_flags(MDArena *arena, int nflags, const char *flags)
+{
+ if (arena == NULL)
+ return 0;
+
+ if (nflags == 0 || flags == NULL) {
+ if (arena->alchem_flags != NULL)
+ free(arena->alchem_flags);
+ arena->alchem_flags = NULL;
+ arena->nalchem_flags = 0;
+ return 0;
+ }
+
+ if (arena->alchem_flags == NULL)
+ arena->alchem_flags = (char *)malloc(nflags);
+ else arena->alchem_flags = (char *)realloc(arena->alchem_flags, nflags);
+ if (arena->alchem_flags == NULL) {
+ arena->nalchem_flags = 0;
+ return -1;
+ }
+ memmove(arena->alchem_flags, flags, nflags);
+ return 0;
+}
+
const char *
md_prepare(MDArena *arena, int check_only)
{
if (t2 > 0)
md_log(arena, "Number of fixed atoms = %d\n", t2);
-
if (arena->natoms_uniq < mol->natoms) {
md_log(arena, "Number of symmetry-unique atoms = %d\n", arena->natoms_uniq);
t2 = 0;
md_log(arena, " %11s %11s %11s %11s %11s %11s %11s %11s", "VDW", "ELECT", "AUX", "SURFACE", "KINETIC", "NET", "TEMP", "TEMP_AVG");
if (periodic)
md_log(arena, " %11s", "VOLUME");
+ if (arena->nalchem_flags > 0)
+ md_log(arena, " %11s %11s", "LAMBDA", "DEDL");
md_log(arena, "\n");
}
md_log(arena, "ENERGY: %11d %11.5f", arena->step, arena->total_energy * INTERNAL2KCAL);
VecCross(v, *av, *bv);
md_log(arena, " %11.5f", VecDot(v, *cv));
}
+ if (arena->nalchem_flags > 0) {
+ md_log(arena, " %11.5f %11.5f", arena->alchem_lambda, arena->alchem_energy / arena->alchem_dlambda);
+ }
md_log(arena, "\n");
md_log(arena, NULL);
}
free(arena->custom_bond_pars);
if (arena->custom_pars != NULL)
free(arena->custom_pars);
+ if (arena->alchem_flags != NULL)
+ free(arena->alchem_flags);
if (arena->bond_par_i != NULL)
free(arena->bond_par_i);
if (arena->angle_par_i != NULL)
Int ncustom_pars;
MDCustomPar *custom_pars;
+ /* Alchemical perturbation */
+ Int nalchem_flags;
+ char *alchem_flags; /* 0: normal, 1: vanishing, 2: appearing */
+ Double alchem_lambda;
+ Double alchem_dlambda;
+ Double alchem_energy;
+
/* Fix atoms */
/* Int nfix_atoms;
fix_atom_record *fix_atoms; */
int md_check_abnormal_dihedral(MDArena *arena, Molecule *mol, int idx);
int md_check_abnormal_improper(MDArena *arena, Molecule *mol, int idx);
+int md_set_alchemical_flags(MDArena *arena, int nflags, const char *flags);
+
void md_amend_by_symmetry(MDArena *arena);
void md_cell_recalculate(MDArena *arena);
void md_scale_cell(MDArena *arena, const Transform tf, int scale_atoms);
int i;
Vector *vdr;
Double limit, elimit, dielec_r;
+ Double lambda, dlambda;
+
MDVerlet *vl;
/* MDVdwCache *vdw_cache = arena->vdw_cache; */
Atom *atoms = arena->mol->atoms;
&& !(get_group_flag(group_flags_1, vl->n2) && get_group_flag(group_flags_2, vl->n1)))
continue;
}
+
+ if (arena->nalchem_flags > 0) {
+ char c1, c2;
+ if (vl->n1 < arena->nalchem_flags)
+ c1 = arena->alchem_flags[vl->n1];
+ else c1 = 0;
+ if (vl->n2 < arena->nalchem_flags)
+ c2 = arena->alchem_flags[vl->n2];
+ else c2 = 0;
+ if ((c1 == 1 && c2 == 2) || (c1 == 2 && c2 == 1))
+ continue;
+ if (c1 == 1 || c2 == 1) {
+ lambda = (1.0 - arena->alchem_lambda);
+ dlambda = -arena->alchem_dlambda;
+ } else if (c1 == 2 || c2 == 2) {
+ lambda = arena->alchem_lambda;
+ dlambda = arena->alchem_dlambda;
+ } else {
+ lambda = 1.0;
+ dlambda = 0.0;
+ }
+ } else {
+ lambda = 1.0;
+ dlambda = 0.0;
+ }
+
if (vl->vdw_type == 1) {
A = vp->par.A14;
B = vp->par.B14;
k0 *= arena->scale14_vdw;
k1 *= arena->scale14_vdw;
}
+ k0 /= vl->mult;
+ k1 *= lambda;
VecScale(fij, rij, k1);
- *energies += k0 / vl->mult;
+ *energies += k0 * lambda;
if (forces != NULL) {
VecDec(forces[vl->n1], fij);
VecInc(forces[vl->n2], fij);
}
+ if (dlambda != 0.0)
+ arena->alchem_energy += k0 * dlambda;
+
if (arena->debug_result && arena->debug_output_level > 1) {
fprintf(arena->debug_result, "nonbonded(vdw) force %d-%d: r=%f, k0=%f, k1=%f, {%f %f %f}\n", vl->n1+1, vl->n2+1, sqrt(r2), k0/KCAL2INTERNAL, k1*sqrt(r2)/KCAL2INTERNAL, fij.x/KCAL2INTERNAL, fij.y/KCAL2INTERNAL, fij.z/KCAL2INTERNAL);
}
memset(arena->energies, 0, sizeof(Double) * kSlowIndex);
memset(arena->forces, 0, sizeof(Vector) * kSlowIndex * natoms);
arena->total_energy = 0.0;
-
+ arena->alchem_energy = 0.0;
+
if (arena->step == arena->start_step || arena->step % arena->surface_potential_freq == 0) {
doSurface = 1;
arena->energies[kSurfaceIndex] = 0.0;
s_AverageTempSym, s_AndersenFreqSym, s_AndersenCouplingSym, s_RandomSeedSym,
s_DielectricSym, s_GradientConvergenceSym, s_CoordinateConvergenceSym, s_UseXplorShiftSym,
s_Scale14VdwSym, s_Scale14ElectSym, s_RelocateCenterSym, s_SurfaceProbeRadiusSym,
-s_SurfaceTensionSym, s_SurfacePotentialFreqSym, s_UseGraphiteSym;
+s_SurfaceTensionSym, s_SurfacePotentialFreqSym, s_UseGraphiteSym,
+s_AlchemicalLambdaSym, s_AlchemicalDeltaLambdaSym;
struct s_MDArenaAttrDef {
char *name;
{"surface_tension", &s_SurfaceTensionSym, 0, 0, 'f', offsetof(MDArena, surface_tension)},
{"surface_potential_freq", &s_SurfacePotentialFreqSym, 0, 0, 'i', offsetof(MDArena, surface_potential_freq)},
{"use_graphite", &s_UseGraphiteSym, 0, 0, 'i', offsetof(MDArena, use_graphite)},
+ {"alchemical_lambda", &s_AlchemicalLambdaSym, 0, 0, 'f', offsetof(MDArena, alchem_lambda)},
+ {"alchemical_delta_lambda", &s_AlchemicalDeltaLambdaSym, 0, 0, 'f', offsetof(MDArena, alchem_dlambda)},
{NULL} /* Sentinel */
};
/*
* call-seq:
+ * set_alchemical_perturbation(group1, group2) -> [group1, group2]
+ *
+ * Set vanishing and appearing atom groups for alchemical perturbation.
+ */
+static VALUE
+s_MDArena_SetAlchemicalPerturbation(VALUE self, VALUE gval1, VALUE gval2)
+{
+ IntGroup *ig1, *ig2;
+ MDArena *arena;
+ char *flags;
+ int i, n;
+ Data_Get_Struct(self, MDArena, arena);
+ if (gval1 == Qnil && gval2 == Qnil) {
+ md_set_alchemical_flags(arena, 0, NULL);
+ return self;
+ }
+ if (arena->mol == NULL)
+ rb_raise(rb_eMolbyError, "Molecule is not set");
+ n = arena->mol->natoms;
+ flags = (char *)calloc(1, n);
+ ig1 = IntGroupFromValue(gval1); /* nil argument is taken as an empty IntGroup */
+ ig2 = IntGroupFromValue(gval2);
+ for (i = 0; i < n; i++) {
+ if (IntGroupLookupPoint(ig1, i) >= 0)
+ flags[i] = 1;
+ if (IntGroupLookupPoint(ig2, i) >= 0) {
+ if (flags[i] == 1)
+ rb_raise(rb_eMolbyError, "duplicate atom (%d) in vanishing and appearing groups", i);
+ flags[i] = 2;
+ }
+ }
+ if (md_set_alchemical_flags(arena, n, flags) != 0)
+ rb_raise(rb_eMolbyError, "cannot set alchemical flags");
+ free(flags);
+ gval1 = ValueFromIntGroup(ig1);
+ gval2 = ValueFromIntGroup(ig2);
+ IntGroupRelease(ig1);
+ IntGroupRelease(ig2);
+ return rb_ary_new3(2, gval1, gval2);
+}
+
+/*
+ * call-seq:
+ * get_alchemical_perturbation -> [group1, group2] or nil
+ *
+ * If alchemical perturbation is enabled, get the current vanishing and appearing atom groups.
+ * Otherwise, return nil.
+ */
+static VALUE
+s_MDArena_GetAlchemicalPerturbation(VALUE self)
+{
+ IntGroup *ig1, *ig2;
+ VALUE gval1, gval2;
+ MDArena *arena;
+ int i;
+ Data_Get_Struct(self, MDArena, arena);
+ if (arena->nalchem_flags == 0)
+ return Qnil;
+ if (arena->mol == NULL)
+ rb_raise(rb_eMolbyError, "Molecule is not set");
+ ig1 = IntGroupNew();
+ ig2 = IntGroupNew();
+ for (i = 0; i < arena->nalchem_flags; i++) {
+ if (arena->alchem_flags[i] == 1)
+ IntGroupAdd(ig1, i, 1);
+ else if (arena->alchem_flags[i] == 2)
+ IntGroupAdd(ig2, i, 1);
+ }
+ gval1 = ValueFromIntGroup(ig1);
+ gval2 = ValueFromIntGroup(ig2);
+ IntGroupRelease(ig1);
+ IntGroupRelease(ig2);
+ return rb_ary_new3(2, gval1, gval2);
+}
+
+/*
+ * call-seq:
* keys -> Array
*
* Returns an array of valid attributes.
rb_define_method(rb_cMDArena, "[]=", s_MDArena_Set, 2);
rb_define_method(rb_cMDArena, "to_hash", s_MDArena_ToHash, 0);
rb_define_method(rb_cMDArena, "keys", s_MDArena_Keys, 0);
+ rb_define_method(rb_cMDArena, "set_alchemical_perturbation", s_MDArena_SetAlchemicalPerturbation, 2);
+ rb_define_method(rb_cMDArena, "get_alchemical_perturbation", s_MDArena_GetAlchemicalPerturbation, 0);
rb_define_method(rb_cMDArena, "print_surface_area", s_MDArena_PrintSurfaceArea, 0);
/* All setter and getter are handled with the same C function (attribute name is taken
OSDN Git Service
