Temperature surges simulating calcite surface with water

Temperature surges simulating calcite surface with water

Dear Lammps users,

I have been trying to simulate a surface of calcium carbonate surface in contact with water. Using the force field published by Raiteri et al in J. Phys Chem C 2015, 119, 24447-24458 DOI:10.1021/acs.jpcc.5b07532
The version of Lammps I am currently using for testing is the Ubuntu binary, obtained in ppa:gladky-anton/lammps
The structure data file was created using Molteplate. The dynamics computation details are the same as published in the above mentioned paper (as accurately as described in the article) The full input files are attached.
To facilitate the reading here are my main input parameters

units metal
atom_style full
boundary p p p
read_data system.data
include forcefield
variable Nrun equal 10000
variable Nf equal ${Nrun}/1000
variable Ne equal 10
variable Nr equal ${Nf}/${Ne}
variable Text equal 300.0
velocity all create 300 1234
velocity all scale 300.00
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
fix integrate all nvt temp ${Text} ${Text} 0.1 tchain 5
compute T all temp
fix TempAve all ave/time ${Ne} ${Nr} ${Nf} c_T
thermo_style custom time temp f_TempAve etotal pe enthalpy cella
thermo_modify flush yes
thermo ${Nf}
dump trj all atom 100 sys.lammpstrj
dump_modify trj scale no sort id
run ${Nrun}

With this, the initial lines of the thermo output are:

Setting up Verlet run ...
Unit style : metal
Current step : 0
Time step : 0.001
Per MPI rank memory allocation (min/avg/max) = 38.53 | 38.53 | 38.53 Mbytes
Time Temp f_TempAve TotEng PotEng Enthalpy Cella
0 300 0 -31404.299 -31733.136 -31917.052 42.186
0.1 1029.5886 795.24135 -27208.552 -28337.111 -28067.012 42.186
0.2 758.85354 832.2868 -25304.993 -26136.793 -28280.761 42.186
0.3 790.00988 817.5759 -24610.583 -25476.533 -22370.756 42.186
0.4 752.93261 762.70952 -24623.035 -25448.344 -23777.292 42.186
0.5 645.94172 696.76808 -24416.257 -25124.291 -26312.212 42.186

I have been searching the mailing list, but finding working input files of surfaces is a bit difficult. I have also scanned the manual, changed the thermostat to Berendsen, that one controls the temperature very well, but similar system was published before with Nosé-Hoover. The paper cited uses npt, I have done so as well with the same results, but I have tried to reduce the problem
I would like some advice regarding the jump in the temperature from 300 K to over 1000K. It will gradually return to 300K but that initial surge causes huge relaxations that are not expected nor seeing with other codes.
I have tried other systems Ca, CO3 ions in solution. The temperature behaves very well in that case, but as a calcite surface is being published before, made me think I am missing something important in my input file.
Please note the run time used here is for testing purposes only.

Best regards,
Sergio

Dear Axel, 

 

Thank you for your quick reply, 

 

I am going to try to reply in line as you did. 

________________________________________
From: Axel Kohlmeyer [akohlmey@gmail.com]
Sent: 26 April 2018 15:29
To: Ruiz Hernandez, S.E. (Sergio)
Cc: lammps-users@lists.sourceforge.net
Subject: Re: [lammps-users] Temperature surges simulating calcite surface with water

On Thu, Apr 26, 2018 at 7:13 AM, Ruiz Hernandez, S.E. (Sergio)
<s.e.ruizhernandez@uu.nl> wrote:
> Dear Lammps users,
>
> I have been trying to simulate a surface of calcium carbonate surface in contact with water. Using the force field published by Raiteri et al in J. Phys Chem C 2015, 119, 24447-24458 DOI:10.1021/acs.jpcc.5b07532
> The version of Lammps I am currently using for testing is the Ubuntu binary, obtained in ppa:gladky-anton/lammps

that is a useless description. as explained in (the beginning of) the
LAMMPS manual, the version of LAMMPS is the date printed at the
beginning of every run.

 

Here it is,  LAMMPS (23 Oct 2017)

 

 

> The structure data file was created using Molteplate. The dynamics computation details are the same as published in the above mentioned paper (as accurately as described in the article) The full input files are attached.

[...]

> I have been searching the mailing list, but finding working input files of surfaces is a bit difficult. I have also scanned the manual, changed the thermostat to Berendsen, that one controls the temperature very well, but similar system was published before with Nosé-Hoover. The paper cited uses npt, I have done so as well with the same results, but I have tried to reduce the problem
> I would like some advice regarding the jump in the temperature from 300 K to over 1000K. It will gradually return to 300K but that initial surge causes huge relaxations that are not expected nor seeing with other codes.

are you referring to running the *exact same* input geometry and force
field with a different MD software?

 

Yes indeed same structure, same force field and DL_POLY 4. Sorry but If I was referring to Microsoft Paint and AMBER FF I would not be asking on this list or bothering you.   

> I have tried other systems Ca, CO3 ions in solution. The temperature behaves very well in that case, but as a calcite surface is being published before, made me think I am missing something important in my input file.

you seem to be missing that a) almost any simulation protocol needs
include some steps of equilibration, and b) that many times initial
geometries may not be fully in accord with the force field used, even
if they represent experimental structures.

 

I am not missing (a) on purpose. but thanks to pointing that out, I cannot find anywhere in the manual instructions that are explicit to equilibration period. Other than start from lower temperatures, NVE or drastic thermostat controls like "drag" I have tried several suggestions from the manual to ease the starting of the simulation without success. And please note these structure is being already simulated in DL_POLY and is behaving as expected. Having said so, (b) was done. I would really appreciate the page of the manual in which equilibrations methods are treated.        

thus for many complex systems, simulation protocols have to start with
checks for overlaps and high potential energy areas, usually followed
minimization/geometry optimization, followed by a
equilibration/thermalization MD runs. for complex systems and delicate
structures, multiple minimizations alternating with short
equilibration runs may be needed, often also for each subsystem
separately.

 

if you have a polar surface in contact with a polar liquid, it is
*extremely* likely, that there are high potential energy spots where
the two subsystems are in contact, since the tool you are using to
construct the system does not no anything about the strength of the
interactions between the two systems.

 

Oh, I see I was a bit misleading with my use of MOLTEMPLATE. The tool was used to translate a pdb geometry of the last geometry of the system that was produced by DL_POLY (same FF)  Sorry 

furthermore, before starting a simulation with a multi-component
system, you should validate your force field parameter and geometry
input by simulating each component separately and verify, that there
are no typos or unit conversion errors or geometry mistakes.

 

at this point i do feel a bit like a broken record, since i seem to be
giving this kind of advice and suggestion regularly on this mailing
list, so perhaps you need to widen your scope of search through the
mailing list archives a bit.

axel.

 

I am really sorry to make you feel like a broken record. As I said before, the force field has been used before in this system. I am using a FF file that was provided in a Supplementary Information of a paper and, I am actually using the values in another well known MD code. In fact, I have simulated the water alone, the ions on their own and both systems together using the same version of lammps. The issues ONLY arise when it is a surface and, the system is not exploding or going totally crazy. It is just relaxing far too much (in comparison) and showing initial spikes of temperature that disappear after the first few steps . It is happening with  Nosé-Hoover thermostat not with the others. 

 

I am sure to be missing some very important parameter to equilibrate the system. I am 100% sure there are not numerical errors in the FF and the initial geometry is not far from equilibrium. I am saying this because according to what is published the calcite surface and the bulk behaved very well in lammps. My appeal is for help of finding the right keyword to fire the simulation.

 

Best wishes,

Sergio 

>
> Best regards,
> Sergio
>