Skip to content

Commit

Permalink
Implemented direct output of CVs and other internal PLUMED variables
Browse files Browse the repository at this point in the history
Use plumed cmd, and extras to accumulate them
  • Loading branch information
ceriottm committed Apr 20, 2024
1 parent 4b8c618 commit f110482
Show file tree
Hide file tree
Showing 14 changed files with 202 additions and 290 deletions.
26 changes: 23 additions & 3 deletions ipi/engine/forcefields.py
Original file line number Diff line number Diff line change
Expand Up @@ -644,7 +644,6 @@ def evaluate(self, r):
r["status"] = "Done"
r["t_finished"] = time.time()


class FFPlumed(FFEval):
"""Direct PLUMED interface
Expand All @@ -671,6 +670,7 @@ def __init__(
init_file="",
plumeddat="",
plumedstep=0,
plumedextras=[],
):
"""Initialises FFPlumed.
Expand All @@ -689,6 +689,7 @@ def __init__(
self.plumed = plumed.Plumed()
self.plumeddat = plumeddat
self.plumedstep = plumedstep
self.plumedextras = plumedextras
self.init_file = init_file

if self.init_file.mode == "xyz":
Expand Down Expand Up @@ -718,6 +719,21 @@ def __init__(
self.plumedrestart = True
self.plumed.cmd("setRestart", 1)
self.plumed.cmd("init")

ipi_name = lambda name : "ipi_"+name.replace(".","_")
self.plumed_data = {}
for x in plumedextras:
rank = np.zeros( 1, dtype=np.int_ )
self.plumed.cmd(f"getDataRank {x}", rank )
if rank[0]>1:
raise ValueError("Cannot retrieve varibles with rank > 1")
shape = np.zeros( rank[0], dtype=np.int_ )
if shape[0]>1:
raise ValueError("Cannot retrieve varibles with size > 1")
self.plumed.cmd(f"getDataShape {x}", shape )
self.plumed_data[x] = np.zeros(shape, dtype=np.double)
self.plumed.cmd(f"setMemoryForData {x}", self.plumed_data[x] )

self.charges = dstrip(myatoms.q) * 0.0
self.masses = dstrip(myatoms.m)
self.lastq = np.zeros(3 * self.natoms)
Expand Down Expand Up @@ -759,7 +775,7 @@ def evaluate(self, r):
self.plumed.cmd("setVirial", vir)
self.plumed.cmd("prepareCalc")
self.plumed.cmd("performCalcNoUpdate")

bias = np.zeros(1, float)
self.plumed.cmd("getBias", bias)
v = bias[0]
Expand All @@ -771,8 +787,12 @@ def evaluate(self, r):
f[:] -= dstrip(self.system_force.f).flatten()
vir[:] -= -dstrip(self.system_force.vir)

extras = {"raw": ""}
for x in self.plumed_data:
extras[str(x)] = self.plumed_data[x].copy()

# nb: the virial is a symmetric tensor, so we don't need to transpose
r["result"] = [v, f, vir, {"raw": ""}]
r["result"] = [v, f, vir, extras]
r["status"] = "Done"

def mtd_update(self, pos, cell):
Expand Down
9 changes: 5 additions & 4 deletions ipi/engine/forces.py
Original file line number Diff line number Diff line change
Expand Up @@ -400,7 +400,7 @@ def bind(self, beads, cell, fflist, output_maker):
)
self._extras = depend_value(
name="extras",
value=np.zeros(self.nbeads, float),
value={},
func=self.extra_gather,
dependencies=[self._forces[b]._extra for b in range(self.nbeads)],
)
Expand Down Expand Up @@ -530,6 +530,7 @@ class ScaledForceComponent:
def __init__(self, baseforce, scaling=1):
self.bf = baseforce
self.name = baseforce.name
self.nbeads = baseforce.nbeads
self.ffield = baseforce.ffield
self._scaling = depend_value(name="scaling", value=scaling)
self._f = depend_array(
Expand All @@ -554,10 +555,10 @@ def __init__(self, baseforce, scaling=1):
value=np.zeros((self.bf.nbeads, 3, 3)),
dependencies=[self.bf._virs, self._scaling],
)
self._extras = depend_array(
self._extras = depend_value(
name="extras",
func=lambda: self.bf.extras,
value=np.zeros(self.bf.nbeads),
value={},
dependencies=[self.bf._extras],
)

Expand Down Expand Up @@ -1447,7 +1448,7 @@ def extra_combine(self):

re = {}
for k in range(self.nforces):
# combines the extras from the different force components
# combines the extras from the different force components
for e, v in self.mforces[k].extras.items():
if e in self.mforces[k].force_extras:
# extras that are tagged as force_extras are treated exactly as if they were an energy/force/stress
Expand Down
5 changes: 3 additions & 2 deletions ipi/engine/outputs.py
Original file line number Diff line number Diff line change
Expand Up @@ -346,13 +346,14 @@ def open_stream(self, mode):
"forces",
"extras",
"extras_component",
"extras_bias",
"forces_sc",
"momenta",
]:
# must write out trajectories for each bead, so must create b streams

# prepare format string for file name
if getkey(self.what) == "extras" or getkey(self.what) == "extras_component":
if getkey(self.what)[:6] == "extras":
fmt_fn = self.filename + "_" + fmt_bead
else:
fmt_fn = self.filename + "_" + fmt_bead + "." + self.format
Expand Down Expand Up @@ -477,7 +478,7 @@ def write_traj(
"""

key = getkey(what)
if key in ["extras", "extras_component"]:
if key in ["extras", "extras_component", "extras_bias"]:
stream.write(
" #%s(%s)# Step: %10d Bead: %5d \n"
% (key.upper(), self.extra_type, self.system.simul.step + 1, b)
Expand Down
19 changes: 4 additions & 15 deletions ipi/engine/properties.py
Original file line number Diff line number Diff line change
Expand Up @@ -1122,21 +1122,6 @@ def get_kincv(self, atom=""):
acv = np.dot(q.flatten(), f.flatten())
acv *= -0.5 / self.beads.nbeads
acv += ncount * 1.5 * Constants.kb * self.ensemble.temp
# ~ acv = 0.0
# ~ ncount = 0
# ~
# ~ for i in range(self.beads.natoms):
# ~ if (atom != "" and iatom != i and latom != self.beads.names[i]):
# ~ continue
# ~
# ~ kcv = 0.0
# ~ k = 3*i
# ~ for b in range(self.beads.nbeads):
# ~ kcv += q[b,k]* f[b,k] + q[b,k+1]* f[b,k+1] + q[b,k+2]* f[b,k+2]
# ~ kcv *= -0.5/self.beads.nbeads
# ~ kcv += 1.5*Constants.kb*self.ensemble.temp
# ~ acv += kcv
# ~ ncount += 1

if ncount == 0:
# TODO: don't warn if bosons are matched
Expand Down Expand Up @@ -2853,6 +2838,10 @@ def __init__(self):
Fetches the extras from a specific force component, indicated in parentheses [extras_component(idx)]. """,
"func": (lambda idx: self.system.forces.extras_component(int(idx))),
},
"extras_bias": {
"help": """The additional data returned by the bias forcefield, printed verbatim or expanded as a dictionary. See "extras". """,
"func": (lambda: self.system.ensemble.bias.extras),
},
"isotope_zetatd": {
"dimension": "undefined",
"help": """Thermodynamic isotope fractionation direct estimator in the form of ratios of partition functions. Takes two arguments, 'alpha' , which gives the
Expand Down
16 changes: 14 additions & 2 deletions ipi/inputs/forcefields.py
Original file line number Diff line number Diff line change
Expand Up @@ -521,14 +521,24 @@ class InputFFPlumed(InputForceField):
"help": "The current step counter for PLUMED calls",
},
),
"plumedextras": (
InputArray,
{"dtype": str, "default": input_default(factory=np.zeros, args=(0, str)),
"help": "List of variables defined in the PLUMED input, that should be transferred to i-PI as `extras` fields." },
),
}

attribs = {}

attribs.update(InputForceField.attribs)
fields.update(InputForceField.fields)

default_help = """ Direct PLUMED interface """
default_help = """
Direct PLUMED interface. Can be used to implement metadynamics in i-PI in combination with
the <metad> SMotion class. NB: if you use PLUMED for constant biasing (e.g. for umbrella sampling)
the bias will be computed but there will be no output as specified in the `plumed.dat` file
unless you include a <metad> tag, that triggers the log update. """

default_label = "FFPLUMED"

def store(self, ff):
Expand All @@ -539,6 +549,7 @@ def store(self, ff):
# self.plumedstep.store(pstep)
self.plumedstep.store(ff.plumedstep)
self.file.store(ff.init_file)
self.plumedextras.store(np.array(ff.plumed_data.keys()))

def fetch(self):
super(InputFFPlumed, self).fetch()
Expand All @@ -551,7 +562,8 @@ def fetch(self):
threaded=self.threaded.fetch(),
plumeddat=self.plumeddat.fetch(),
plumedstep=self.plumedstep.fetch(),
init_file=self.file.fetch(),
plumedextras=self.plumedextras.fetch(),
init_file=self.file.fetch(),
)


Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -3,4 +3,4 @@
ref_simulation.frc_c.xyz xyz
ref_simulation.pos_c.xyz xyz
ref_simulation.out numpy
ref_COLVAR numpy
ref_simulation.cv_0 numpy
8 changes: 4 additions & 4 deletions ipi_tests/regression_tests/tests/PLUMED/BIAS.NOAUTO/input.xml
Original file line number Diff line number Diff line change
@@ -1,8 +1,10 @@
<simulation mode="static" verbosity="medium">
<output prefix='simulation'>
<properties stride='1' filename='out'> [ step, potential{electronvolt}] </properties>
<properties stride='10' filename='out'> [ step, potential{electronvolt}, ensemble_bias{electronvolt}] </properties>
<trajectory stride="20" filename="pos_c" format="xyz"> x_centroid </trajectory>
<trajectory stride="20" filename="frc_c" format="xyz"> f_centroid </trajectory>
<trajectory stride="10" filename="cv" extra_type="dhh"> extras_bias </trajectory>
<trajectory stride="10" filename="rest" extra_type="rest.bias"> extras_bias </trajectory>
</output>
<total_steps> 100 </total_steps>
<ffsocket name="pswater" mode="unix">
Expand All @@ -11,6 +13,7 @@
<ffplumed name="plumed">
<file mode="xyz">./init.xyz</file>
<plumeddat> plumed.dat </plumeddat>
<plumedextras> [dhh, rest.bias] </plumedextras>
</ffplumed>
<system >
<initialize nbeads='1'>
Expand All @@ -35,7 +38,4 @@
</dynamics>
</motion>
</system>
<smotion mode="metad">
<metad> <metaff> [ plumed ] </metaff> </metad>
</smotion>
</simulation>
26 changes: 0 additions & 26 deletions ipi_tests/regression_tests/tests/PLUMED/BIAS.NOAUTO/ref_COLVAR

This file was deleted.

113 changes: 12 additions & 101 deletions ipi_tests/regression_tests/tests/PLUMED/BIAS.NOAUTO/ref_simulation.out
Original file line number Diff line number Diff line change
@@ -1,103 +1,14 @@
# column 1 --> step : The current simulation time step.
# column 2 --> potential{electronvolt} : The physical system potential energy.
0.00000000e+00 3.25615492e-02
1.00000000e+00 5.76949314e-02
2.00000000e+00 1.41165251e-01
3.00000000e+00 3.66274711e-01
4.00000000e+00 7.39957788e-01
5.00000000e+00 1.18657709e+00
6.00000000e+00 1.64059978e+00
7.00000000e+00 2.11730348e+00
8.00000000e+00 2.51283801e+00
9.00000000e+00 2.80217377e+00
1.00000000e+01 2.98085379e+00
1.10000000e+01 3.02318306e+00
1.20000000e+01 2.90767615e+00
1.30000000e+01 2.71609007e+00
1.40000000e+01 2.45544674e+00
1.50000000e+01 2.13952149e+00
1.60000000e+01 1.78209900e+00
1.70000000e+01 1.45127851e+00
1.80000000e+01 1.18472483e+00
1.90000000e+01 9.72369496e-01
2.00000000e+01 8.25580153e-01
2.10000000e+01 7.68242936e-01
2.20000000e+01 7.84240488e-01
2.30000000e+01 9.09298701e-01
2.40000000e+01 1.13361847e+00
2.50000000e+01 1.38850907e+00
2.60000000e+01 1.64658535e+00
2.70000000e+01 1.88173226e+00
2.80000000e+01 2.05674122e+00
2.90000000e+01 2.19192921e+00
3.00000000e+01 2.28081628e+00
3.10000000e+01 2.34212600e+00
3.20000000e+01 2.34394085e+00
3.30000000e+01 2.30719281e+00
3.40000000e+01 2.13131605e+00
3.50000000e+01 1.96317480e+00
3.60000000e+01 1.74668285e+00
3.70000000e+01 1.50102140e+00
3.80000000e+01 1.30438679e+00
3.90000000e+01 1.12391496e+00
4.00000000e+01 9.69291798e-01
4.10000000e+01 8.66893892e-01
4.20000000e+01 7.85790720e-01
4.30000000e+01 7.61510117e-01
4.40000000e+01 7.99245224e-01
4.50000000e+01 8.98409910e-01
4.60000000e+01 1.04454649e+00
4.70000000e+01 1.22633875e+00
4.80000000e+01 1.40156356e+00
4.90000000e+01 1.56699046e+00
5.00000000e+01 1.74743959e+00
5.10000000e+01 1.89622664e+00
5.20000000e+01 1.99519129e+00
5.30000000e+01 2.04454885e+00
5.40000000e+01 2.05902285e+00
5.50000000e+01 2.05184435e+00
5.60000000e+01 2.03091989e+00
5.70000000e+01 1.96961817e+00
5.80000000e+01 1.89824694e+00
5.90000000e+01 1.78203375e+00
6.00000000e+01 1.62576266e+00
6.10000000e+01 1.43878440e+00
6.20000000e+01 1.27078130e+00
6.30000000e+01 1.12019035e+00
6.40000000e+01 1.05533955e+00
6.50000000e+01 1.08185462e+00
6.60000000e+01 1.14155728e+00
6.70000000e+01 1.20894642e+00
6.80000000e+01 1.25352795e+00
6.90000000e+01 1.33414251e+00
7.00000000e+01 1.39149390e+00
7.10000000e+01 1.47005568e+00
7.20000000e+01 1.53239691e+00
7.30000000e+01 1.58728429e+00
7.40000000e+01 1.59932992e+00
7.50000000e+01 1.60508486e+00
7.60000000e+01 1.59916236e+00
7.70000000e+01 1.55469308e+00
7.80000000e+01 1.52514794e+00
7.90000000e+01 1.49101237e+00
8.00000000e+01 1.44091236e+00
8.10000000e+01 1.41637975e+00
8.20000000e+01 1.39748798e+00
8.30000000e+01 1.34828554e+00
8.40000000e+01 1.28372097e+00
8.50000000e+01 1.24049201e+00
8.60000000e+01 1.17313816e+00
8.70000000e+01 1.12352079e+00
8.80000000e+01 1.11917209e+00
8.90000000e+01 1.10723717e+00
9.00000000e+01 1.14348946e+00
9.10000000e+01 1.17896254e+00
9.20000000e+01 1.25654189e+00
9.30000000e+01 1.34171864e+00
9.40000000e+01 1.43036501e+00
9.50000000e+01 1.49167593e+00
9.60000000e+01 1.51978189e+00
9.70000000e+01 1.51103238e+00
9.80000000e+01 1.46516604e+00
9.90000000e+01 1.40203871e+00
1.00000000e+02 1.37586058e+00
# column 3 --> ensemble_bias{electronvolt} : The bias applied to the current ensemble
0.00000000e+00 3.25615492e-02 4.31033085e+00
1.00000000e+01 2.98085379e+00 2.30964185e-01
2.00000000e+01 8.25580153e-01 2.10170212e+00
3.00000000e+01 2.28081628e+00 5.14931223e-01
4.00000000e+01 9.69291798e-01 1.63948920e+00
5.00000000e+01 1.74743959e+00 9.00071457e-01
6.00000000e+01 1.62576266e+00 9.73399614e-01
7.00000000e+01 1.39149390e+00 1.15067255e+00
8.00000000e+01 1.44091236e+00 1.10196105e+00
9.00000000e+01 1.14348946e+00 1.41818509e+00
1.00000000e+02 1.37586058e+00 1.16128765e+00
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
----------------------------------------------------------
ref_simulation.frc_c.xyz xyz
ref_simulation.pos_c.xyz xyz
ref_simulation.dhh_0 numpy
ref_simulation.out numpy
ref_COLVAR numpy
ref_HILLS numpy
Loading

0 comments on commit f110482

Please sign in to comment.