search

peteboyd / lammps_interface

automatic generation of LAMMPS input files for molecular dynamics simulations of MOFs
MIT License
125 stars 63 forks source link

Inconsistent atom types from published work and from lammps_interface.py #15

Closed askforarun closed 5 years ago

askforarun commented 5 years ago

Hi Peter.

Thank you for sharing the code. I tried to redo some of the simulations published in your recent article Force-Field Prediction of Materials Properties in Metal-OrganicFrameworks

I started with IRMOF-1 series. I converted the data file data. IRMOF-1.repeat to cif using topotools , obabel. Then I used the lammps_interface.py to obtain the data file.

I found the atom types corresponding to Dreiding forcefield used in the original paper does not match the one obtained by using lammps_interface.py.

Published work (IRMOF-1. repeat Direiding) Step E_vdwl E_coul E_pair E_bond E_angle E_dihed E_impro 0 506.73555 -14322.243 -13815.507 678.80669 209.03348 2.0872193e-13 15360

From lammps_interface.py

Step E_vdwl E_coul E_pair E_bond E_angle E_dihed E_impro 0 506.80088 0 506.80088 708.09748 280.49672 1.9914626e-13 0

I am more interested in the inconsistency between angle, dihedral, improper energies. one additional angle type is missing in the second one and O3 is identified as O2.

It would be helpful if you could clarify this Arun

askforarun commented 5 years ago

cif from ata.IRMOF-1.repeat from published work. I used topo tools to read the data file, converted to pdb and then from pdb to cif using mercury. IRMOF-1.docx If you use lammps_interface.py on this (with dreiding forcefield) it should return the same data file as the published one.

However atom O3 is identified as O2 and the angle, dihedral energies dont match as shown.

Thank you again Arun

peteboyd commented 5 years ago

Hi Arun, There are several reasons that the code will not produce the same answer as the paper produced two years ago. From experience, force fields are rather fragile and must be used with the understanding that the results can change with different (yet equally valid) interpretations. For UFF and Dreiding, they were not designed with coordination polymers in mind, thus atom typing was subject to some interpretation by the authors (particularly where metal coordination geometry was concerned).

Through the steps you took to obtain a .cif file from the data file, it is possible the code has interpreted the atom types differently. I'm not entirely sure what topotools and obabel will do to the bonding structure and atom types, but it is possible that this step has caused the changes in atom typing you observe. In addition, the code has undergone numerous changes over the years, with users noticing bugs and suggesting fixes that will ultimately change the interpretation of a cif file. You can check the commit history of this code to get an understanding of all the changes.

For these reasons I am not surprised you are getting a different interpretation of the force field. What would be more concerning is if the results of a calculation are strange or unphysical. In which case you should send me all of the information you have used to generate the lammps input files and tell me the version of Lammps you are using to run the calculation. One of the reasons we provide the data files with the article was so that, at the time when the article was published, one could reproduce the results as they were calculated then.

Pete

askforarun commented 5 years ago

Thank you for your reply Peter

Arun