ab initio REMD exchanges

[Martín Taccone]

Hi everyone! I'm running some ab initio REMD with i-PI and FHI-aims of a DNA-Ag system of ~130 atoms to try to explore some conformational space. Due to computational resources I've ran a 10 ps NVT dynamic with 4 replicas (from 300K to 900K) for testing. Although the best choice would obviously be to increase replicas and simulation time, it probably be too much. But my question arise from the fact that when I follow the total energy from FHI-aims output as a function of time I get what it seems to look like a lot of exchanges between replicas, but when I follow i-PI total energy of the ensembles (and swap file for exchanges) I only see a few exchanges between replicas. Am I missing something? Why there is such difference between both cases? I leave my i-Pi input and two images of ensembles total energy.

<simulation verbosity='low'>
  <total_time>  3.02400e+07</total_time>
  <total_steps>10000</total_steps>
  <ffsocket name='fhiaims'>
    <latency> 1.0000000e-02 </latency>
    <slots>4</slots>
    <timeout> 8.00000000e+04</timeout>
    <address>ipihost</address>
    <port> 14123 </port>
  </ffsocket>
  <output prefix='test'>
    <properties stride='1' filename='out'>  [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}] </properties>
    <trajectory filename='pos' stride='1' format='xyz'> positions{angstrom} </trajectory>
    <trajectory filename='vel' stride='1' format='xyz'> velocities </trajectory>
    <trajectory filename='xc' stride='1' format='pdb'> x_centroid{angstrom} </trajectory>
    <trajectory filename='for' stride='1' format='xyz'> forces </trajectory>
    <checkpoint stride='200' filename='chk' overwrite="false"/>
    <checkpoint stride='2' filename='restart' overwrite="true"/>
  </output>
  <prng>
    <seed>31415</seed>
  </prng>
  <system prefix='300K'>
    <forces>
      <force forcefield='fhiaims'> </force>
    </forces>
    <initialize nbeads='1'>
      <file mode='xyz'> init.xyz </file>
      <cell mode='abc'> [200.0, 200.0, 200.0] </cell>
      <velocities mode='thermal' units='kelvin'> 300 </velocities>
    </initialize>
    <ensemble>
      <temperature units='kelvin'> 300 </temperature>
    </ensemble>
    <motion mode='dynamics'>
      <dynamics mode='nvt'>
        <timestep units='femtosecond'> 1.0 </timestep>
    <thermostat mode="langevin">
     <tau units="femtosecond"> 50 </tau>
    </thermostat>
      </dynamics>
    </motion>
  </system>
  <system prefix='500K'>
    <forces>
      <force forcefield='fhiaims'> </force>
    </forces>
    <initialize nbeads='1'>
      <file mode='xyz'> init.xyz </file>
      <cell mode='abc'> [200.0, 200.0, 200.0] </cell>
      <velocities mode='thermal' units='kelvin'> 500 </velocities>
    </initialize>
    <ensemble>
      <temperature units='kelvin'> 500 </temperature>
    </ensemble>
    <motion mode='dynamics'>
      <dynamics mode='nvt'>
        <timestep units='femtosecond'> 1.0 </timestep>
    <thermostat mode="langevin">
     <tau units="femtosecond"> 50 </tau>
    </thermostat>
      </dynamics>
    </motion>
  </system>
  <system prefix='700K'>
    <forces>
      <force forcefield='fhiaims'> </force>
    </forces>
    <initialize nbeads='1'>
      <file mode='xyz'> init.xyz </file>
      <cell mode='abc'> [200.0, 200.0, 200.0] </cell>
      <velocities mode='thermal' units='kelvin'> 700 </velocities>
    </initialize>
    <ensemble>
      <temperature units='kelvin'> 700 </temperature>
    </ensemble>
    <motion mode='dynamics'>
      <dynamics mode='nvt'>
        <timestep units='femtosecond'> 1.0 </timestep>
    <thermostat mode="langevin">
     <tau units="femtosecond"> 50 </tau>
    </thermostat>
      </dynamics>
    </motion>
  </system>
  <system prefix='900K'>
    <forces>
      <force forcefield='fhiaims'> </force>
    </forces>
    <initialize nbeads='1'>
      <file mode='xyz'> init.xyz </file>
      <cell mode='abc'> [200.0, 200.0, 200.0] </cell>
      <velocities mode='thermal' units='kelvin'> 900 </velocities>
    </initialize>
    <ensemble>
      <temperature units='kelvin'> 900 </temperature>
    </ensemble>
    <motion mode='dynamics'>
      <dynamics mode='nvt'>
        <timestep units='femtosecond'> 1.0 </timestep>
    <thermostat mode="langevin">
     <tau units="femtosecond"> 50 </tau>
    </thermostat>
      </dynamics>
    </motion>
  </system>
  <smotion mode="remd">
   <remd>
    <stride> 10 </stride>
    <swapfile> exchanges.out </swapfile>
   </remd>
  </smotion>
</simulation>

[Michele Ceriotti]

Is the trajectory "re-ordered"? After finishing the trajectory you should run the post-processing step 

i-pi-remdsort input.xml

that will re-print the trajectories so that each trajectory correspond to one ensemble. 

For many reasons, the implementation of REMD in i-PI exchanges the ensembles, not

the configurations, so each output file transitions between different ensembles. It is usually a good 

idea to print in the output ensemble_temp, that will tell you what is the ensemble at each point

in time without having to re-order the trajectory.

Hope this clarifies your problem.

Best

Michele

[Martín Taccone]

The trayectory was not re-ordered. Ok, that makes everything more clear.

Thanks!

Martín.