Is it possible to perform NPT simulation with schnetpack?

[siddharth-ncl-work]

Hi, I have gone through the schnetpack tutorials but couldn't find NPT simulation example (let me know if I missed it). Is it possible to perform NPT simulation with schnetpack? If yes, can you please explain how to do so. Thanks

[mgastegger]

It is possible to run NPT simulations with the master branch, although this is not documented yet.

To do so, you have to specify a barostat in the dynamics section of the input, for example:

dynamics:
  n_steps: 100000
  integrator:
    type: verlet 
    time_step: 0.50
  barostat:
    type: nhc_iso
    target_pressure: 1.0 # bar
    temperature: 300 # K
    time_constant: 100 # fs
    time_constant_cell: 1000 # fs
    time_constant_barostat: 500 # fs
   massive: true

This would use a Martyna-Tobias-Klein barostat for isotropic cell fluctuations. In this case, no additional thermostat needs to be specified, since particle and cell temperatures are controlled via NHC thermostats.

In addition to the barostat, the calculator used for the MD also needs to compute the stress tensors, e.g. with the SchNetPack calculator, this can be done via:

calculator:
  type: schnet
  model_file: <PATH/TO/MODEL>
  required_properties:
  - energy
  - forces
  - stress  # < for stress
  cutoff: <MODEL CUTOFF>
  force_handle: forces
  stress_handle: stress # < for stress
  position_conversion: Angstrom
  force_conversion: eV/Angstrom
  stress_conversion: eV/A/A/A # < for stress
  neighbor_list: ase

Finally, a simulation cell and PBC need to be provided in the structure input, this means the comment line of the extended XYZ file would read something like this:

Lattice="34.982601165771484 0.0 0.0 0.0 34.982601165771484 0.0 0.0 0.0 4.982601165771484" Properties=species:S:1:pos:R:3 pbc="T T T"

[siddharth-ncl-work]

@mgastegger thanks a lot for getting back to me. I will give it a try on ethanol and update you on the results.

I have a few doubts regarding simulation setup:

1. For NPT, do we also need a model trained on stress tensors along with the energy and forces? I followed this tutorial and have a trained model for ethanol that gives energy and forces....will that work?

2. I am working on molecular system of benzene and methanol (will be trying on ethanol as well), I think simulation cell (ie. box) would be required but do I also need to specify PBC (ie. lattice vectors)?

3. Is massive: true part of barostat: or calculator: block? I think there is some issue with its indentation.

4. How to monitor pressure and volume of the system?

[siddharth-ncl-work]

[UPDATE] I was able to fix the following error by reinstalling SchNetPack from the cloned repo. I don't need help with this error but I would appreciate if you could clarify the doubts from earlier comment

I tried running NPT simulation on ethanol as per the instructions but got following error. I used spk_md.py file for the simulation. The repo was cloned again to make sure I was on the latest commit of the master branch. SchNet model was obtained by following this tutorial.

../../../src/scripts/spk_md.py:14: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
  tradoffs = yaml.load(tf)
INFO:root:Read options from input_script_NPT.
INFO:root:Random state initialized with seed 2752731239
INFO:root:Create model directory
/home/vanka/anaconda3/envs/torchmd/lib/python3.8/site-packages/ase/atoms.py:967: VisibleDeprecationWarning: Use get_global_number_of_atoms() instead
  warnings.warn('Use get_global_number_of_atoms() instead',
INFO:root:Loaded model from best_model
Traceback (most recent call last):
  File "../../../src/scripts/spk_md.py", line 27, in <module>
    md = MDSimulation(config)
  File "/home/vanka/anaconda3/envs/torchmd/lib/python3.8/site-packages/schnetpack/md/parsers/md_setup.py", line 59, in __init__
    SetupCalculator(self)
  File "/home/vanka/anaconda3/envs/torchmd/lib/python3.8/site-packages/schnetpack/md/parsers/md_setup.py", line 127, in __init__
    self._setup(md_initializer)
  File "/home/vanka/anaconda3/envs/torchmd/lib/python3.8/site-packages/schnetpack/md/parsers/md_setup.py", line 355, in _setup
    calculator = CalculatorInit(calculator_dict).initialized
  File "/home/vanka/anaconda3/envs/torchmd/lib/python3.8/site-packages/schnetpack/md/parsers/md_options.py", line 181, in __init__
    self.initialized = target_class(*inputs, **extra_kwargs)
TypeError: __init__() got an unexpected keyword argument 'stress_handle'

Input xyz file (ethanol_with_pbc.xyz)

9
Lattice="10.0 0.0 0.0 0.0 10.0 0.0 0.0 0.0 10.0" Properties=species:S:1:pos:R:3 pbc="T T T"
O        3.36190000       5.44035000       5.00300000
C        4.48680000       4.57415000       5.00300000
C        5.75060000       5.40745000       5.00300000
H        4.43730000       3.92865000       5.88490000
H        4.43790000       3.94685000       4.10840000
H        6.63810000       4.76865000       4.98530000
H        5.77570000       6.07135000       4.13260000
H        5.79470000       6.04585000       5.89160000
H        3.40400000       5.97705000       5.81290000

Input Script

device: cuda
simulation_dir: test_output
overwrite: True
dynamics:
  n_steps: 100000
  integrator:
    type: verlet
    time_step: 0.50
  barostat:
    type: nhc_iso
    target_pressure: 1.0 # bar
    temperature: 300 # K
    time_constant: 100 # fs
    time_constant_cell: 1000 # fs
    time_constant_barostat: 500 # fs
    massive: true
calculator:
  type: schnet
  model_file: best_model #model from tutorial
  required_properties:
  - energy
  - forces
  - stress # < for stress
  force_handle: forces
  stress_handle: stress # < for stress
  position_conversion: Angstrom
  force_conversion: kcal/mol/Angstrom
  stress_conversion: kcal/mol/A/A/A
system:
  molecule_file: ethanol_with_pbc.xyz
  n_replicas: 1
  initializer:
    type: maxwell-boltzmann
    temperature: 300
    remove_translation: true
    remove_rotation: true
logging:
  file_logger:
    buffer_size: 100
    streams:
    - molecules
    - properties
  write_checkpoints: 100

Could you please have look at this?

Thank you

[mgastegger]

Once again sorry for the late reply. I am glad to hear you could solve the stress_handle error.

Regarding your questions:

1. We are currently testing this ourselves. So far it seems that including the energies and forces is enough to also get sufficiently accurate stress tensor predictions.
2. The periodic boundaries need to be specified in the extended XYZ file as well. For example pbc="T T T means that PBCs are applied in x, y and z directions.
3. There was an indentation missing. massive: true should be part of the barostat block.
4. In the current schnetpack version, is should be possible to monitor the pressure by adding pressure_logger: 1 to the logging block. Unfortunately, the volume needs to be extracted from the HDF5 file. We are, however. currently finalizing a new version of SchNetPack with better support for NPT simulations.