Cg minimization + shake algorithm

LAMMPS LAMMPS Beginners
1

Dear LAMMPS users
I am trying to implement a lammps code from a paper about boiling of water on cooper plate. I,ve written the below code but I have a problem. In the paper min_style cg and SHAKE algorithm are used. When I use both of them the simulation process after minimization does NOT proceed, and when I comment min_style cg , the simulation diverge. The main part of the code is as bellow. All data are based on a paper and are correct.
I appreciate it if anyone help me with this issue.

Thanks
Hamed

#--------------------------------------------------------------#

Intalization

#--------------------------------------------------------------#
dimension 3
boundary p p p
units real
atom_style full

#--------------------------------------------------------------#

Atom Definition

#--------------------------------------------------------------#
read_data data.txt
region Cu block 0 122.91 0 21.69 0 29.00 units box
lattice fcc 3.615
create_atoms 3 region Cu

group H type 1
group O type 2
group Water type 1 2
group Cu type 3

set type 1 charge 0.52
set type 2 charge -1.04
set type 3 charge 0.00

#--------------------------------------------------------------#

Interatomic Potential

#--------------------------------------------------------------#
pair_style hybrid lj/cut/tip4p/long 2 1 1 1 0.1546 15.0 lj/cut 15
kspace_style pppm/tip4p 1.0e-5

pair_coeff 1 1 lj/cut/tip4p/long 0.0 0.0
pair_coeff 1 2 lj/cut/tip4p/long 0.0 0.0
pair_coeff 2 2 lj/cut/tip4p/long 0.15535 3.166
pair_coeff 1 3 lj/cut 0.0 0.0
pair_coeff 2 3 lj/cut 0.8563 2.54785
pair_coeff 3 3 lj/cut 4.72 1.9297

bond_style harmonic
bond_coeff 1 0.0 0.9572
angle_style harmonic
angle_coeff 1 0.0 104.52

#--------------------------------------------------------------#

Setting and Output of Simulation

#--------------------------------------------------------------#
variable Ne equal 10
variable Nr equal 10
variable Nf equal 100

variable Tstart equal 370.0

neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes

min_style cg
minimize 1.0e-01 1.0e-01 10000 10000

fix CuSpring Cu spring/self 2.9332
fix constrain Water shake 1.0e-4 100 0 b 1 a 1
fix removeMomentum Water momentum 1 linear 1 1 1

fix 7 Water nvt temp {Tstart} {Tstart} 100
fix 8 Water langevin {Tstart} {Tstart} 100 654843

timestep 1.0

thermo_style custom step temp
thermo 1
dump 1 all xyz 1 dump.xyz

run 500000

2

Look at the paper very carefully. Did they use LAMMPS? Did they explicitly state that they were using fix shake during minimization (or only during MD). The SHAKE algorithm is only derived for MD simulation time stepping and cannot be applied during minimization (applying constraints during minimization is a “hard” problem). All MD codes I know that would allow to define bond/angle/etc. constraints during minimization did so by turning off the SHAKE algorithm and instead using a harmonic potential with a very large force constant. LAMMPS is honest in this case and does not trick you. It just refuses to run. It is easy enough to implement the workaround, though.

That is a very bold statement and likely to be proven wrong. After many years of responding to statements similar to yours, I would be willing to wager a sizable amount of money that this statement is incorrect, i.e. that either you are not following the paper exactly or that the information in the paper is inaccurate (or plain wrong). In my personal experience most people (about 95%) that had problems yet claimed their input was perfect, their input turned out to be at fault after all.

3

Thanks for your time and responese Axel,

From the paper: “All the simulations were performed using the open source classical MD code Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)”

From the paper: “The hydrogen‑oxygen bonds and angle of water molecules are kept rigid using the SHAKE algorithm.”
and
“At the beginning of each simulation, the energy of the entire system was minimized using the conjugate gradient (CG) algorithm.”

4

That is rather vague.

5

So how can I first minimize using cg, and then run my simulation with shake algorithm?
Is it possible at all?
I’m literally confused :innocent:

6

I already told you how.

7

I’m sorry for asking this again…I’m almost new to LAMMPS simulation, would you please do me a favor and explain this to me more clearer?

Thanks

8

What is confusing about using a large force constant?

9

My question is how can I first do minimization while the Shake algorithm is turned off? and then how can I turn on the shake to run my simulation?

10

This should be obvious from studying the LAMMPS manual and specifically the sections explaining how input files and the individual commands (e.g. fix) work. If you want to be using LAMMPS in a meaningful way, you have to understand how its commands work and how it processes input. Do not expect people like me give you “for X do this, not that” kind of instructions.

11

Thank you Axel for your time and consideration.

Hamed