Fix units (for good)

bin/i-pi-mux-positions

@@ -5,6 +5,7 @@ from __future__ import print_function
 import os
 import sys
 import argparse
+from ipi.utils.messages import verbosity
 
 from ipi.utils.io import open_backup, iter_file_name, print_file
 
@@ -19,6 +20,8 @@ input files.
 
 def main(fns_in, fn_out, begin, end, stride):
 
+    verbosity.level = "low" 
+
     print('Multiplexing {:d} beads into one trajectory.'.format(len(fns_in)))
     print()
 

doc/manual.tex

@@ -470,6 +470,15 @@ \subsection{Internal units}
 }
 
 
+Regarding the specification of these units in the i-PI input files, the user is able to specify units both in the i-PI input file or in the structure file. 
+If the structure file is of the .xyz format, the units specifications should be present in the comment line. Examples of such inputs can be found in \texttt{examples/pes-regtest/io-units/} .
+The code then behaves in the following way, depending on the user's choice: 
+\begin{itemize}
+\item If no units are specified in the input, i-PI tries to guess from the comment line of the structure file and it nothing is present, assumes atomic units, or angstrom for PDB files.
+\item If units are specified in both input and structure file and they match, conversion happens just once. If they do not match, an error is raised and i-PI stops.
+\end{itemize}
+
+
 \section{Core features}
 
 The functionality of \ipi includes:

drivers/driver.f90

@@ -222,7 +222,7 @@ PROGRAM DRIVER
       ELSEIF (vstyle == 11) THEN
          IF (par_count .ne. 3) THEN
             WRITE(*,*) "Error:  incorrect initialization string included for qtip4pf-efield. &
-          Provide the three components of the electric field in V/nm"
+     &    Provide the three components of the electric field in V/nm"
             STOP "ENDED"
       ELSE
             ! We take in an electric field in volts / nm.This must be converted 

examples/lammps/h2o-piglet.8/input.xml

@@ -2,6 +2,8 @@
     <output prefix="simulation">
         <properties stride="1" filename="out">  [ step, time{picosecond}, conserved{kelvin}, temperature{kelvin}, kinetic_cv{kelvin}, potential{kelvin}, pressure_cv{megapascal} ] </properties>
         <trajectory filename="pos" stride="20"> positions </trajectory>
+        <trajectory filename="kin" stride="20"> kinetic_cv </trajectory>
+        <trajectory filename="kod" stride="20"> kinetic_od </trajectory>
     </output>
     <total_steps> 500000 </total_steps>
     <prng><seed> 32342 </seed></prng>

examples/lammps/h2o-pimd+mts.4/input.xml

@@ -1,7 +1,7 @@
 <simulation mode="md">
     <output prefix="simulation">
         <properties stride="1" filename="out">  [ step, time{picosecond}, conserved, temperature{kelvin}, kinetic_cv, potential, pressure_cv{megapascal} ] </properties>
-        <trajectory filename="pos" stride="20"> positions </trajectory>
+        <trajectory filename="pos" stride="1"> positions </trajectory>
         <checkpoint stride="200"/>
     </output>
     <total_steps>10000</total_steps>

examples/pes-regtest/io-units/h2o-angstrom.xyz

@@ -0,0 +1,5 @@
+3
+# CELL(abcABC):   18.89726    18.89726    18.89726    90.00000    90.00000    90.00000  Step:          40  Bead:       2 positions{angstrom}  cell{atomic_unit}
+       O  4.97336e+00  4.92337e+00  4.96921e+00
+       H  5.71284e+00  5.57553e+00  4.90940e+00
+       H  5.01842e+00  4.53605e+00  4.05255e+00

examples/pes-regtest/io-units/h2o-nocomment.xyz

@@ -0,0 +1,5 @@
+3
+
+O          5.00000        5.00000        5.00000
+H          5.98343        5.00000        5.00000
+H          4.84583        4.376419       4.10217

examples/pes-regtest/io-units/h2o-velocities.xyz

@@ -0,0 +1,5 @@
+3
+# CELL(abcABC):   18.89726    18.89726    18.89726    90.00000    90.00000    90.00000  Step:          40  Bead:       2 velocities{atomic_unit}  cell{atomic_unit}
+       O -1.57474e-04 -6.43328e-04 -3.02692e-04
+       H  4.75131e-04  1.92827e-03 -9.58797e-04
+       H  1.42102e-03  2.09636e-03  2.30678e-03

examples/pes-regtest/io-units/input.xml

@@ -0,0 +1,46 @@
+<!-- REGTEST
+DEPENDENCIES  h2o-nocomment.xyz h2o-velocities.xyz h2o-angstrom.xyz
+COMMAND    i-pi-driver -h localhost -p 31400 -m qtip4pf
+ENDREGTEST -->
+
+<simulation verbosity="high">
+  <output prefix='PSmono'>
+    <properties stride='1' filename='out'>  [ time{picosecond}, conserved, temperature{kelvin}, kinetic_md, kinetic_cv, kinetic_cv(O), kinetic_cv(H), potential ] </properties>
+    <trajectory filename='pos' format="xyz" stride='2'> positions{angstrom} </trajectory>
+    <trajectory filename='vel' format="xyz" stride='5'> velocities </trajectory>
+    <trajectory filename='for' format="xyz" stride='2'> forces </trajectory>
+    <trajectory filename='kin' format="xyz" stride='2'> kinetic_cv </trajectory>
+    <trajectory filename='dip' stride='2'> extras </trajectory>
+    <checkpoint stride='100' filename='chk' overwrite='true'/>
+    <checkpoint stride='2000' filename='restart' overwrite='false'/>
+  </output>
+  <total_steps>40</total_steps>
+  <prng><seed>23658</seed></prng>
+  <ffsocket mode='inet' pbc='false' name='driver'>
+    <address>localhost</address>
+    <!-- <address>driver</address> -->
+      <port>31400</port> <latency>0.02</latency> <timeout>400</timeout>
+  </ffsocket>
+  <system>
+    <initialize nbeads='4'>
+      <file mode='xyz' units='angstrom'> h2o-nocomment.xyz </file>
+      <cell mode="abc" units="angstrom"> [ 10.0, 10.0, 10.0 ] </cell>
+      <velocities mode="thermal" units="kelvin"> 300 </velocities>
+      <velocities mode='xyz' bead='3'> h2o-velocities.xyz </velocities>
+      <positions mode='xyz' bead='2'> h2o-angstrom.xyz </positions>
+    </initialize>
+    <forces><force forcefield='driver'/></forces>
+    <ensemble>
+      <temperature units="kelvin">300</temperature>
+    </ensemble>
+    <motion mode='dynamics'>
+      <fixcom>True</fixcom>
+      <dynamics mode='nvt'>
+	<timestep units="femtosecond">0.5</timestep>
+	<thermostat mode='langevin'>
+	  <tau units='femtosecond'> 100</tau>
+	</thermostat>
+      </dynamics>
+    </motion>
+  </system>
+</simulation>

examples/pes-regtest/zundel/input.xml

@@ -1,8 +1,8 @@
-<!--REGTEST
+<!--REGTEST 
 DEPENDENCIES  h5o2.dms4B.coeff.com.dat h5o2.pes4B.coeff.dat h5o2+.xyz
 COMMAND(10)    i-pi-driver -u -h REGTEST_SOCKET -m zundel
 COMMAND(6)    i-pi-driver -u -h REGTEST_SOCKET -m zundel
-ENDREGTEST-->
+ ENDREGTEST-->
 
 <simulation>
     <ffsocket mode='unix' name='driver'>

ipi/engine/forces.py

@@ -731,7 +731,6 @@ def forces_4th_order(self, index):
          warning("ERROR: Suzuki-Chin factorization requires even number of beads!")
          exit()
 
-
       # calculates the finite displacement.
       fbase = depstrip(self.f)
       eps = self.mforces[index].epsilon
@@ -950,4 +949,3 @@ def get_coeffsc_part_2(self):
       rc[0::2] =  (self.alpha / self.omegan2 / 9.0)
       rc[1::2] =  ((1.0 - self.alpha) / self.omegan2 / 9.0)
       return np.asmatrix(rc).T
-

ipi/engine/initializer.py

@@ -26,7 +26,7 @@
 from ipi.utils.messages import verbosity, warning, info
 
 
-__all__ = ['Initializer', 'InitBase', 'InitIndexed']
+__all__ = ['Initializer', 'InitBase', 'InitIndexed', 'InitFile']
 
 
 class InitBase(dobject):
@@ -70,8 +70,8 @@ class InitIndexed(InitBase):
       bead: Which bead to initialize the value of.
    """
 
-   def __init__(self, value="", mode="", units="", index=-1, bead=-1):
-      """Initializes InitFile.
+   def __init__(self, value="", mode="", units="", index=-1, bead=-1, **others):
+      """Initializes InitIndexed.
 
       Args:
          value: A string which specifies what value to initialize the
@@ -85,8 +85,23 @@ def __init__(self, value="", mode="", units="", index=-1, bead=-1):
 
       super(InitIndexed,self).__init__(value=value, mode=mode, units=units, index=index, bead=bead)
 
-
-def init_file(mode, filename):
+class InitFile(InitBase):
+    def __init__(self, value="", mode="", units="", cell_units="", bead=-1, **others):
+        """Initializes InitIndexed.
+
+        Args:
+            value: A string which specifies what value to initialize the
+                simulation property to.
+            mode: A string specifiying what style of initialization should be
+                used to read the data.
+            units: A string giving which unit the value is in.
+            cell_units: A string giving which unit the cell parameters for the files are
+            bead: Which bead to initialize the value of.
+        """
+
+        super(InitFile,self).__init__(value=value, mode=mode, units=units, cell_units=cell_units, bead=bead)
+
+def init_file(mode, filename, dimension="length", units="automatic", cell_units="automatic"):
    """Reads a @mode file and returns the data contained in it.
 
    Args:
@@ -100,11 +115,13 @@ def init_file(mode, filename):
 
    rfile = open(filename, "r")
    ratoms = []
+
+   info("Initializing from file %s. Dimension: %s, units: %s, cell_units: %s" % (filename, dimension, units, cell_units), verbosity.low)
    while True:
    #while loop, so that more than one configuration can be given
    #so multiple beads can be initialized at once.
       try:
-         ret = read_file(mode, rfile)
+         ret = read_file(mode, rfile, dimension=dimension, units=units, cell_units=cell_units)
       except EOFError:
          break
       ratoms.append(ret["atoms"])
@@ -139,7 +156,7 @@ def init_chk(filename):
    return (rbeads, rcell, rmotion)
 
 
-def init_beads(iif, nbeads):
+def init_beads(iif, nbeads, dimension="length", units="automatic", cell_units="automatic"):
    """Initializes a beads object from an appropriate initializer object.
 
    Args:
@@ -157,7 +174,7 @@ def init_beads(iif, nbeads):
    elif mode == "manual":
       raise ValueError("Cannot initialize manually a whole beads object.")
    else:
-      ret = init_file(mode, value)
+      ret = init_file(mode, value, dimension, units, cell_units)
       ratoms = ret[0]
       rbeads = Beads(ratoms[0].natoms,len(ratoms))
       for i in range(len(ratoms)):
@@ -166,7 +183,7 @@ def init_beads(iif, nbeads):
    return rbeads
 
 
-def init_vector(iif, nbeads, momenta=False):
+def init_vector(iif, nbeads, momenta=False, dimension="length", units="automatic", cell_units="automatic"):
    """Initializes a vector from an appropriate initializer object.
 
    Args:
@@ -179,7 +196,7 @@ def init_vector(iif, nbeads, momenta=False):
    mode = iif.mode
    value = iif.value
    if mode == "xyz" or mode == "pdb":
-      rq = init_beads(iif, nbeads).q
+      rq = init_beads(iif, nbeads, dimension, units, cell_units).q
    elif mode == "chk":
       if momenta: rq = init_beads(iif, nbeads).p
       else:       rq = init_beads(iif, nbeads).q
@@ -310,7 +327,7 @@ def init_stage1(self, simul):
 
          if k == "cell":
             if v.mode == "manual":
-                rh = v.value.reshape((3,3))
+                rh = v.value.reshape((3,3))*unit_to_internal("length",v.units,1.0)
             elif v.mode == "chk":
                rh = init_chk(v.value)[1].h
             elif init_file(v.mode,v.value)[1].h.trace() == -3:
@@ -319,15 +336,15 @@ def init_stage1(self, simul):
                #+set to -1 from the io_units and nothing is read here.
                continue
             else:
-               rh =  init_file(v.mode,v.value)[1].h
+               rh =  init_file(v.mode,v.value,cell_units=v.units)[1].h
 
             if fcell :
                warning("Overwriting previous cell parameters", verbosity.low)
 
 
-            warning_units_message(v, 'cell')
+            #warning_units_message(v, 'cell')
 
-            rh *= unit_to_internal("length",v.units,1.0)
+            #rh *= unit_to_internal("length",v.units,1.0)
 
             simul.cell.h = rh
             if simul.cell.V == 0.0:
@@ -340,10 +357,9 @@ def init_stage1(self, simul):
             if fmass:
                warning("Overwriting previous atomic masses", verbosity.medium)
             if v.mode == "manual":
-               rm = v.value
+               rm = v.value * unit_to_internal("mass",v.units,1.0)
             else:
                rm = init_beads(v, self.nbeads).m
-            rm *= unit_to_internal("mass",v.units,1.0)
 
             if v.bead < 0:   # we are initializing the path
                if (fmom and fmass):
@@ -382,11 +398,9 @@ def init_stage1(self, simul):
             if fpos:
                warning("Overwriting previous atomic positions", verbosity.medium)
             # read the atomic positions as a vector
-            rq = init_vector(v, self.nbeads)
-
-            warning_units_message(v, 'position')
+
+            rq = init_vector(v, self.nbeads, dimension="length", units=v.units)
 
-            rq *= unit_to_internal("length",v.units,1.0)
             nbeads, natoms = rq.shape
             natoms /= 3
 
@@ -446,8 +460,7 @@ def init_stage1(self, simul):
             if fmom:
                warning("Overwriting previous atomic momenta", verbosity.medium)
             # read the atomic momenta as a vector
-            rp = init_vector(v, self.nbeads, momenta=True)
-            rp *= unit_to_internal("momentum", v.units, 1.0)
+            rp = init_vector(v, self.nbeads, momenta=True, dimension="momentum", units=v.units)
             nbeads, natoms = rp.shape
             natoms /= 3
 
@@ -461,14 +474,13 @@ def init_stage1(self, simul):
             set_vector(v, simul.beads.p, rp)
             fmom = True
 
-            warning_units_message(v, 'momenta')
+            #warning_units_message(v, 'momenta')
 
          elif k == "velocities":
             if fmom:
                warning("Overwriting previous atomic momenta", verbosity.medium)
             # read the atomic velocities as a vector
-            rv = init_vector(v, self.nbeads)
-            rv *= unit_to_internal("velocity", v.units, 1.0)
+            rv = init_vector(v, self.nbeads, dimension="velocity", units=v.units)
             nbeads, natoms = rv.shape
             natoms /= 3
 
@@ -487,7 +499,7 @@ def init_stage1(self, simul):
             set_vector(v, simul.beads.p, rv)
             fmom = True
 
-            warning_units_message(v, 'velocity')
+            #warning_units_message(v, 'velocity')
 
          elif k == "gle": pass   # thermostats must be initialised in a second stage
 

ipi/engine/motion/dynamics.py

@@ -409,8 +409,8 @@ def step(self, step=None):
         self.ptime += time.time()
 
         self.ttime -= time.time()
-        self.barostat.thermostat.step()
         self.thermostat.step()
+        self.barostat.thermostat.step()
         self.pconstraints()
         self.ttime += time.time()