Pressure on Perovskite cell


I have been trying to simulate pressure on a 2d perovskite. Creating a thermostat was doable by using temp/berendsen. I then wanted to get a constant pressure and afterwards look at the rdf. However I got stuck on the barostat. Both the thermo press output as the f_3 from the press/berendsen did not output the ‘requested’ pressure.
I know that the run times are low at the moment. I also didn’t see them converge at run 10.000 and since I am for now running it on a laptop, I haven’t been comfortable to run it for longer. In a week I will be able to upload it to server with parallel computing powers, then ofcourse I will switch to longer calculating times.

I put all the input files in the following onedrive folder:

units       	real
boundary 		p p p
atom_style  	full

pair_style      hybrid lj/charmm/coul/long 15.2 15.3 buck/coul/long 15.2 15.3
bond_style      harmonic
angle_style     harmonic
dihedral_style  fourier

kspace_style 	pppm 1.0e-4
kspace_modify 	diff ad

read_data   	lmp_atoms
replicate       2 2 1  	

include		paircoeff_new

timestep    	0.1

velocity 	all create 300.0 4540 mom yes dist gaussian

thermo_style    custom step time enthalpy etotal press temp vol
thermo          500

fix             1 all nve
fix             2 all temp/berendsen 300.0 300.0 100.0

run 7000

unfix 2

fix 	3 all press/berendsen iso 1 1 $(100.0*dt)  

thermo_style    custom step time enthalpy etotal press temp vol f_3

run 7000

compute         rad all rdf 500 1 2
fix     4 all ave/time 500 10 5000 c_rad[*] file tmp.rdf mode vector

run 7000

Have you checked the individual pressure tensor components? Specifically Pxx, Pyy, and Pzz? Are they all the same? Total pressure is the average of those three and you are enforcing an isotropic cell deformation, which means that you may be forcing different dimensions into the opposite direction of where they would “want” to go.

Good question, I ran the simulation again. When I tried running again I got a non-numeric box error:

ERROR on proc 0: Non-numeric box dimensions - simulation unstable (src/KSPACE/pppm.cpp:1862)
Last command: run 10000

When replacing press/berendsen and nve for nph. nph seemed to do the same, pressure was still in the thousands, looking at pxx, pyy and pzz they were still higher then the 1 atmosphere I tried to get.
What still baffles me is that I use real units, so my pressure should be outputted in atmosphere (atm). However in the output all pressure is in the thousands, which I think would be really high at just above room temperature (300K).

You have a rather hard material that is not very compressible and thus pressure will change a lot with just small changes of the volume.

There are a few things that you can do about this:

  • Check that your force field settings are good and your initial geometry is proper and has no overlaps or close contacts. Keep in mind that you may see a lot of error/force cancellation due to symmetry that may cause issues as soon as the symmetry is broken due to assigning kinetic energy and starting the MD.
  • Run a minimization with fix box/relax before starting the MD. Even if your geometry is “perfect” (e.g. from X-ray data) it may not be consistent with your force field. A minimization may relax your system to a lower potential energy and reduce the problems when starting (far) away from the minimum. To break symmetries you may use the displace_atoms command to apply small random displacements.
  • For fix press/berendsen you may need to change the assumed bulk modulus of your system. The default is assuming that of water, IIRC. This is essentially a scaling constant for how much the volume get changed in response to the pressure.
  • For both the Berendsen and the Nose-Hoover barostats, you have to be very careful in choosing the relaxation constant. Due to the reduced compressibility, you must not do large changes to the volume as they will cause very large fluctuations in pressure.
  • You may need to increase the system size if you have a very small system. Pressure fluctuations are (much) larger for smaller systems and thus make it more difficult to get a stable simulation. Ideally you want to get into a slightly overdamped regime where the volume does not fluctuate due to the change of pressure but slowly (but not too slowly) and moves toward a converged target. That is not always easy to do. In the other extreme, you may get resonances and then the simulation will easily become unstable.
  • Since you have long-range Coulomb, you may want to change the kspace convergence parameter to something like 1.0e-6. The 1.0e-4 value is good enough for converged forces, but does not converge pressure that well.
  • If your initial pressure is too high, you may just empirically increase the volume a little bit to lower it.
  • If your compressibility is not isotropic, you will likely need to change to anisotropic box updates.