return -1;
}
memmove(arena->alchem_flags, flags, nflags);
+ arena->nalchem_flags = nflags;
return 0;
}
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, " %11.5f %11.5f", arena->alchem_lambda, arena->alchem_energy * INTERNAL2KCAL / arena->alchem_dlambda);
}
md_log(arena, "\n");
md_log(arena, NULL);
arena->andersen_thermo_freq = 50;
arena->andersen_thermo_coupling = 0.1;
arena->pressure_freq = 0;
+
+ arena->alchem_dlambda = 0.1;
+
/* arena->pressure_coupling = 0.4;
arena->pressure_trial_width = 0.01;
arena->pressure_control_algorithm = 0;
continue; */ /* Non-unique angle */
if (ap[angles[0]].occupancy == 0.0 || ap[angles[1]].occupancy == 0.0 || ap[angles[2]].occupancy == 0.0)
continue; /* Skip non-occupied atoms */
+ if (arena->nalchem_flags > 0) {
+ if (angles[0] < arena->nalchem_flags && angles[1] < arena->nalchem_flags
+ && angles[2] < arena->nalchem_flags
+ && (arena->alchem_flags[angles[0]] | arena->alchem_flags[angles[1]] | arena->alchem_flags[angles[2]]) == 3)
+ 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);
continue; *//* Non-unique dihedral */
if (ap[dihedrals[0]].occupancy == 0.0 || ap[dihedrals[1]].occupancy == 0.0 || ap[dihedrals[2]].occupancy == 0.0 || ap[dihedrals[3]].occupancy == 0.0)
continue; /* Skip non-occupied atoms */
- else tp = dihedral_pars + *dihedral_par_i;
+ if (arena->nalchem_flags > 0) {
+ if (dihedrals[0] < arena->nalchem_flags && dihedrals[1] < arena->nalchem_flags
+ && dihedrals[2] < arena->nalchem_flags && dihedrals[3] < arena->nalchem_flags
+ && (arena->alchem_flags[dihedrals[0]] | arena->alchem_flags[dihedrals[1]]
+ | arena->alchem_flags[dihedrals[2]] | arena->alchem_flags[dihedrals[3]]) == 3)
+ continue; /* Interaction between vanishing and appearing groups is ignored */
+ }
+ tp = dihedral_pars + *dihedral_par_i;
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);
for (i = 0; i < arena->natoms_uniq; i++) {
Vector f;
Double en;
+ Double lambda, dlambda;
Vector v = mol->atoms[i].r;
if (mol->atoms[i].fix_force < 0)
continue;
if (mol->atoms[i].occupancy == 0.0)
continue;
+ if (arena->nalchem_flags > 0) {
+ char c1 = arena->alchem_flags[i];
+ if (c1 == 1) {
+ lambda = (1.0 - arena->alchem_lambda);
+ dlambda = -arena->alchem_dlambda;
+ } else if (c1 == 2) {
+ lambda = arena->alchem_lambda;
+ dlambda = arena->alchem_dlambda;
+ } else {
+ lambda = 1.0;
+ dlambda = 0.0;
+ }
+ } else {
+ lambda = 1.0;
+ dlambda = 0.0;
+ }
en = f.x = f.y = f.z = 0.0;
for (ix = 0; ix < px; ix++) {
for (iy = 0; iy < py; iy++) {
}
}
}
- en *= pxpypz_1;
+ en *= pxpypz_1 * lambda;
VecScaleSelf(f, pxpypz_1);
- graphite->energies[i] = en;
- graphite->last_energy += en;
+ graphite->energies[i] = en * lambda;
+ graphite->last_energy += en * lambda;
graphite->last_forces[i] = f;
+ if (dlambda != 0.0)
+ arena->alchem_energy += en * dlambda;
}
for (i = arena->natoms_uniq; i < mol->natoms; i++) {
Symop symop = mol->atoms[i].symop;
fn = 0;
nframes = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
- if (strstr(buf, "!:atoms") == buf) {
+ if (strncmp(buf, "!:", 2) != 0)
+ continue; /* Skip until section header is found */
+ bufp = buf;
+ strsep(&bufp, " \t\n");
+ if (strcmp(buf, "!:atoms") == 0) {
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
continue;
ap->intCharge = ibuf[3];
}
continue;
- } else if (strstr(buf, "!:atoms_symop") == buf) {
+ } else if (strcmp(buf, "!:atoms_symop") == 0) {
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
i++;
}
continue;
- } else if (strstr(buf, "!:atoms_fix") == buf) {
+ } else if (strcmp(buf, "!:atoms_fix") == 0) {
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
i++;
}
continue;
- } else if (strstr(buf, "!:positions") == buf) {
+ } else if (strcmp(buf, "!:positions") == 0) {
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
mp->nframes = mp->cframe = 0;
}
continue;
- } else if (strstr(buf, "!:bonds") == buf) {
+ } else if (strcmp(buf, "!:bonds") == 0) {
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
continue;
}
}
continue;
- } else if (strstr(buf, "!:angles") == buf) {
+ } else if (strcmp(buf, "!:angles") == 0) {
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
continue;
}
}
continue;
- } else if (strstr(buf, "!:dihedrals") == buf) {
+ } else if (strcmp(buf, "!:dihedrals") == 0) {
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
continue;
}
}
continue;
- } else if (strstr(buf, "!:impropers") == buf) {
+ } else if (strcmp(buf, "!:impropers") == 0) {
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
continue;
}
}
continue;
- } else if (strstr(buf, "!:xtalcell") == buf && mp->cell == NULL) {
+ } else if (strcmp(buf, "!:xtalcell") == 0 && mp->cell == NULL) {
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
continue;
MoleculeSetCell(mp, dbuf[0], dbuf[1], dbuf[2], dbuf[3], dbuf[4], dbuf[5], 0);
}
continue;
- } else if (strstr(buf, "!:symmetry_operations") == buf) {
+ } else if (strcmp(buf, "!:symmetry_operations") == 0) {
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
Transform tr;
}
}
continue;
- } else if (strstr(buf, "!:periodic_box") == buf) {
+ } else if (strcmp(buf, "!:periodic_box") == 0) {
Vector vs[5];
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
}
}
continue;
- } else if (strstr(buf, "!:md_parameters") == buf) {
+ } else if (strcmp(buf, "!:md_parameters") == 0) {
MDArena *arena;
if (mp->arena == NULL)
mp->arena = md_arena_new(NULL);
bufp++;
valp = strsep(&bufp, "\n");
} else valp = NULL;
+ if (strcmp(comp, "alchem_flags") == 0) {
+ j = (valp == NULL ? 0 : atoi(valp));
+ if (j > 0) {
+ valp = (char *)malloc(j);
+ i = 0;
+ while ((k = fgetc(fp)) >= 0) {
+ ungetc(k, fp);
+ if (k < '0' || k > '9') {
+ snprintf(errbuf, errbufsize, "line %d: too few flags in alchem_flags block", lineNumber + 1);
+ free(valp);
+ goto exit;
+ }
+ ReadLine(buf, sizeof buf, fp, &lineNumber);
+ bufp = buf;
+ while (*bufp != 0) {
+ if (*bufp >= '0' && *bufp <= '2') {
+ if (i >= j) {
+ snprintf(errbuf, errbufsize, "line %d: too many flags in alchem_flags block", lineNumber);
+ free(valp);
+ goto exit;
+ }
+ valp[i++] = *bufp - '0';
+ } else if (*bufp != ' ' && *bufp != '\t' && *bufp != '\n') {
+ snprintf(errbuf, errbufsize, "line %d: strange character (0x%02x) in alchem_flags block", lineNumber, (int)*bufp);
+ free(valp);
+ goto exit;
+ }
+ bufp++;
+ }
+ if (i == j)
+ break;
+ }
+ md_set_alchemical_flags(arena, j, valp);
+ free(valp);
+ }
+ continue;
+ }
/* In the following, the redundant "!= NULL" is to suppress suprious warning */
if ((strcmp(comp, "log_file") == 0 && (pp = &arena->log_result_name) != NULL)
|| (strcmp(comp, "coord_file") == 0 && (pp = &arena->coord_result_name) != NULL)
|| (strcmp(comp, "scale14_vdw") == 0 && (dp = &arena->scale14_vdw) != NULL)
|| (strcmp(comp, "scale14_elect") == 0 && (dp = &arena->scale14_elect) != NULL)
|| (strcmp(comp, "surface_probe_radius") == 0 && (dp = &arena->probe_radius) != NULL)
- || (strcmp(comp, "surface_tension") == 0 && (dp = &arena->surface_tension) != NULL)) {
+ || (strcmp(comp, "surface_tension") == 0 && (dp = &arena->surface_tension) != NULL)
+ || (strcmp(comp, "alchemical_lambda") == 0 && (dp = &arena->alchem_lambda) != NULL)
+ || (strcmp(comp, "alchemical_delta_lambda") == 0 && (dp = &arena->alchem_dlambda) != NULL)) {
*dp = (valp == NULL ? 0.0 : strtod(valp, NULL));
}
}
continue;
- } else if (strstr(buf, "!:pressure_control_parameters") == buf) {
+ } else if (strcmp(buf, "!:pressure_control_parameters") == 0) {
MDPressureArena *pressure;
if (mp->arena == NULL)
mp->arena = md_arena_new(mp);
}
}
continue;
- } else if (strstr(buf, "!:velocity") == buf) {
+ } else if (strcmp(buf, "!:velocity") == 0) {
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
i++;
}
continue;
- } else if (strstr(buf, "!:force") == buf) {
+ } else if (strcmp(buf, "!:force") == 0) {
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
i++;
}
continue;
- } else if (strstr(buf, "!:parameter") == buf) {
+ } else if (strcmp(buf, "!:parameter") == 0 || strcmp(buf, "!:parameters") == 0) {
Parameter *par = mp->par;
if (par == NULL) {
mp->par = ParameterNew();
free(bufp);
}
continue;
- } else if (strstr(buf, "!:trackball") == buf) {
+ } else if (strcmp(buf, "!:trackball") == 0) {
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
i++;
}
continue;
- } else if (strstr(buf, "!:view") == buf) {
+ } else if (strcmp(buf, "!:view") == 0) {
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
continue;
fprintf(fp, "surface_tension %g\n", arena->surface_tension);
fprintf(fp, "surface_potential_freq %d\n", arena->surface_potential_freq);
fprintf(fp, "use_graphite %d\n", arena->use_graphite);
+ fprintf(fp, "alchemical_lambda %g\n", arena->alchem_lambda);
+ fprintf(fp, "alchemical_delta_lambda %g\n", arena->alchem_dlambda);
+ if (arena->nalchem_flags > 0) {
+ fprintf(fp, "alchem_flags %d", arena->nalchem_flags);
+ for (i = 0; i < arena->nalchem_flags; i++) {
+ if (i % 60 == 0)
+ fputc('\n', fp);
+ else if (i % 10 == 0)
+ fputc(' ', fp);
+ fputc('0' + arena->alchem_flags[i], fp);
+ }
+ fputc('\n', fp);
+ }
if (arena->pressure != NULL) {
Double *dp;
fprintf(fp, "pressure_freq %d\n", arena->pressure->freq);
}
free(tempv);
mp->nframes = n_new;
- MoleculeSelectFrame(mp, last_inserted, 0);
+ MoleculeSelectFrame(mp, last_inserted, 1);
MoleculeIncrementModifyCount(mp);
__MoleculeUnlock(mp);
return count;
s_DielectricSym, s_GradientConvergenceSym, s_CoordinateConvergenceSym, s_UseXplorShiftSym,
s_Scale14VdwSym, s_Scale14ElectSym, s_RelocateCenterSym, s_SurfaceProbeRadiusSym,
s_SurfaceTensionSym, s_SurfacePotentialFreqSym, s_UseGraphiteSym,
-s_AlchemicalLambdaSym, s_AlchemicalDeltaLambdaSym;
+s_AlchemicalLambdaSym, s_AlchemicalDeltaLambdaSym, s_AlchemicalEnergySym;
struct s_MDArenaAttrDef {
char *name;
VALUE *symref; /* Address of s_LogFileSym etc. */
ID id; /* Ruby ID of the symbol; will be set within Init_MolbyMDTypes() */
ID sid; /* Ruby ID of the symbol plus '='; will be set within Init_MolbyMDTypes() */
- char type; /* s: string (const char *), i: Int, f: Double. Uppercase: read-only. */
+ char type; /* s: string (const char *), i: Int, f: Double, e: Double in energy dimension (unit conversion is necessary). */
+ /* Uppercase: read-only. */
int offset; /* Offset in the MDArena structure. */
};
static struct s_MDArenaAttrDef s_MDArenaAttrDefTable[] = {
{"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)},
+ {"alchemical_energy", &s_AlchemicalEnergySym, 0, 0, 'E', offsetof(MDArena, alchem_energy)},
{NULL} /* Sentinel */
};
case 'f':
case 'F':
return rb_float_new(*((Double *)p));
+ case 'e':
+ case 'E':
+ return rb_float_new(*((Double *)p) * INTERNAL2KCAL);
default:
rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
}
case 'f':
case 'F':
return rb_float_new(*((Double *)pp));
+ case 'e':
+ case 'E':
+ return rb_float_new(*((Double *)pp) * INTERNAL2KCAL);
case 'X':
/* Isotropic pressure only */
return rb_ary_new3(3, rb_float_new(pres->apply[0]), rb_float_new(pres->apply[4]), rb_float_new(pres->apply[8]));
case 'f':
*((Double *)p) = NUM2DBL(rb_Float(val));
return val;
- case 'S': case 'I': case 'F':
+ case 'e':
+ *((Double *)p) = NUM2DBL(rb_Float(val) * KCAL2INTERNAL);
+ return val;
+ case 'S': case 'I': case 'F': case 'E':
rb_raise(rb_eMolbyError, "The attribute '%s' is read-only", rb_id2name(aid));
default:
rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
case 'f':
*((Double *)pp) = NUM2DBL(rb_Float(val));
return val;
+ case 'e':
+ *((Double *)pp) = NUM2DBL(rb_Float(val) * KCAL2INTERNAL);
+ return val;
case 'X':
/* Isotropic pressure only */
val = rb_ary_to_ary(val);
else pres->cell_flexibility[j] = 0.0;
}
return val;
- case 'S': case 'I': case 'F':
+ case 'S': case 'I': case 'F': case 'E':
rb_raise(rb_eMolbyError, "The attribute '%s' is read-only", rb_id2name(aid));
default:
rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
Data_Get_Struct(self, MDArena, arena);
if (gval1 == Qnil && gval2 == Qnil) {
md_set_alchemical_flags(arena, 0, NULL);
- return self;
+ return Qnil;
}
if (arena->mol == NULL)
rb_raise(rb_eMolbyError, "Molecule is not set");
- n = arena->mol->natoms;
+ n = arena->xmol->natoms;
flags = (char *)calloc(1, n);
- ig1 = IntGroupFromValue(gval1); /* nil argument is taken as an empty IntGroup */
- ig2 = IntGroupFromValue(gval2);
+ ig1 = (gval1 == Qnil ? NULL : IntGroupFromValue(gval1));
+ ig2 = (gval2 == Qnil ? NULL : IntGroupFromValue(gval2));
for (i = 0; i < n; i++) {
- if (IntGroupLookupPoint(ig1, i) >= 0)
+ if (ig1 != NULL && IntGroupLookupPoint(ig1, i) >= 0)
flags[i] = 1;
- if (IntGroupLookupPoint(ig2, i) >= 0) {
+ if (ig2 != NULL && 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);
+ if (ig1 != NULL) {
+ gval1 = ValueFromIntGroup(ig1);
+ IntGroupRelease(ig1);
+ }
+ if (ig2 != NULL) {
+ gval2 = ValueFromIntGroup(ig2);
+ IntGroupRelease(ig2);
+ }
return rb_ary_new3(2, gval1, gval2);
}
arena = self.md_arena
keys = arena.to_hash.keys
# The read-only keys
- read_only = [:step, :coord_frame, :transient_temperature, :average_temperature]
+ read_only = [:step, :coord_frame, :transient_temperature, :average_temperature, :alchemical_energy]
# Sort the keys so that they are (a little more) readable
[:timestep, :temperature, :cutoff, :electro_cutoff, :pairlist_distance,
:scale14_vdw, :scale14_elect, :use_xplor_shift, :dielectric,
OSDN Git Service
