[lammps-users] tip4p+shake

majido · February 28, 2022, 9:33pm

Dear lammps users,

I use tip4p and fix shake for my water molecules in desalination but they resist piston pressure (350-500 bar). Water molecules shrink from their original state and behave abnormally. The input script is given below:

dimension 3
boundary p p p
units metal
atom_style full
timestep 0.00005

Atom Definition

read_data membrane.data

group fluid1 id 801:9880
group fluid2 id 10387:11286
group fluid union fluid1 fluid2

group Cg type 1
group Cl type 2
group H type 3
group N type 4
group Na type 5
group O type 6
group S type 7
group Hw type 8
group Ow type 9
group Cm type 10

group H2O union Hw Ow
group salt union Cl Na
group membrane union H N O S Cm

variable zpiston1 equal -130
variable zmembrane equal 0
variable zpiston2 equal 130

variable zmin1 equal {zpiston1}-1.0 variable zmax1 equal {zpiston1}+1.0
region piston1zone block INF INF INF INF {zmin1} {zmax1} units box

variable zmin3 equal {zpiston2}-1.0 variable zmax3 equal {zpiston2}+1.0
region piston2zone block INF INF INF INF {zmin3} {zmax3} units box

group piston1 region piston1zone
group piston2 region piston2zone
group piston union piston1 piston2

group thermostat_target union H2O salt membrane piston1 piston2

set group Cl charge -1
set group Na charge 1
set group Hw charge 0.5242
set group Ow charge -1.0484
set group Cg charge 0

Interatomic Potential

pair_style hybrid lj/cut/tip4p/long 9 8 7 16 0.1546 12 12 airebo 3.0 0 1
pair_modify tail yes
pair_coeff * * airebo CH.airebo C NULL NULL NULL NULL NULL NULL NULL NULL NULL
kspace_style pppm/tip4p 1.0e-4

pair_coeff 2 2 lj/cut/tip4p/long 0.01322 3.48
pair_coeff 3 3 lj/cut/tip4p/long 0.00043 2.85
pair_coeff 4 4 lj/cut/tip4p/long 0.01799 3.56
pair_coeff 5 5 lj/cut/tip4p/long 0.02168 2.80
pair_coeff 6 6 lj/cut/tip4p/long 0.01799 3.30
pair_coeff 7 7 lj/cut/tip4p/long 0.01322 3.77
pair_coeff 8 8 lj/cut/tip4p/long 0.00000 0.0
pair_coeff 9 9 lj/cut/tip4p/long 0.00803 2.81
pair_coeff 10 10 lj/cut/tip4p/long 0.01322 3.72

pair_coeff 1 2 lj/cut/tip4p/long 0.01322 3.60
pair_coeff 1 3 lj/cut/tip4p/long 0.00239 3.28
pair_coeff 1 4 lj/cut/tip4p/long 0.01542 3.64
pair_coeff 1 5 lj/cut/tip4p/long 0.01693 3.26
pair_coeff 1 6 lj/cut/tip4p/long 0.01542 3.51
pair_coeff 1 7 lj/cut/tip4p/long 0.01322 3.75
pair_coeff 1 8 lj/cut/tip4p/long 0.00000 1.86
pair_coeff 1 9 lj/cut/tip4p/long 0.01030 3.27
pair_coeff 1 10 lj/cut/tip4p/long 0.01322 3.72

pair_coeff 2 3 lj/cut/tip4p/long 0.00239 3.17
pair_coeff 2 4 lj/cut/tip4p/long 0.01542 3.52
pair_coeff 2 5 lj/cut/tip4p/long 0.01693 3.14
pair_coeff 2 6 lj/cut/tip4p/long 0.01542 3.39
pair_coeff 2 7 lj/cut/tip4p/long 0.01322 3.63
pair_coeff 2 8 lj/cut/tip4p/long 0.00000 1.74
pair_coeff 2 9 lj/cut/tip4p/long 0.01030 3.15
pair_coeff 2 10 lj/cut/tip4p/long 0.01322 3.60

pair_coeff 3 4 lj/cut/tip4p/long 0.00279 3.20
pair_coeff 3 5 lj/cut/tip4p/long 0.00306 2.82
pair_coeff 3 6 lj/cut/tip4p/long 0.00279 3.07
pair_coeff 3 7 lj/cut/tip4p/long 0.00239 3.31
pair_coeff 3 8 lj/cut/tip4p/long 0.00000 1.42
pair_coeff 3 9 lj/cut/tip4p/long 0.00186 2.83
pair_coeff 3 10 lj/cut/tip4p/long 0.00239 3.28

pair_coeff 4 5 lj/cut/tip4p/long 0.01975 3.18
pair_coeff 4 6 lj/cut/tip4p/long 0.01799 3.43
pair_coeff 4 7 lj/cut/tip4p/long 0.01542 3.66
pair_coeff 4 8 lj/cut/tip4p/long 0.00000 1.78
pair_coeff 4 9 lj/cut/tip4p/long 0.01202 3.18
pair_coeff 4 10 lj/cut/tip4p/long 0.01542 3.64

pair_coeff 5 6 lj/cut/tip4p/long 0.01319 3.05
pair_coeff 5 7 lj/cut/tip4p/long 0.01693 3.29
pair_coeff 5 8 lj/cut/tip4p/long 0.00000 3.26
pair_coeff 5 9 lj/cut/tip4p/long 0.01319 2.80
pair_coeff 5 10 lj/cut/tip4p/long 0.01693 3.26

pair_coeff 6 7 lj/cut/tip4p/long 0.01542 3.54
pair_coeff 6 8 lj/cut/tip4p/long 0.00000 1.65
pair_coeff 6 9 lj/cut/tip4p/long 0.01202 3.06
pair_coeff 6 10 lj/cut/tip4p/long 0.01542 3.51

pair_coeff 7 8 lj/cut/tip4p/long 0.00000 1.88
pair_coeff 7 9 lj/cut/tip4p/long 0.01030 3.29
pair_coeff 7 10 lj/cut/tip4p/long 0.01322 3.75

pair_coeff 8 9 lj/cut/tip4p/long 0.00000 1.40
pair_coeff 8 10 lj/cut/tip4p/long 0.00000 1.86

pair_coeff 9 10 lj/cut/tip4p/long 0.01030 3.27

bond_style harmonic
bond_coeff 1 15.177 1.530
bond_coeff 2 15.177 1.090
bond_coeff 3 15.177 1.462
bond_coeff 4 15.177 1.420
bond_coeff 5 15.177 1.800
bond_coeff 6 15.177 1.022
bond_coeff 7 15.177 0.980
bond_coeff 8 15.177 1.394

angle_style harmonic
angle_coeff 1 2.168 109.471
angle_coeff 2 2.168 109.471
angle_coeff 3 2.168 109.471
angle_coeff 4 2.168 109.471
angle_coeff 5 2.168 109.471
angle_coeff 6 2.168 106.700
angle_coeff 7 2.168 106.700
angle_coeff 8 2.168 104.510
angle_coeff 9 2.168 92.1000
angle_coeff 10 2.168 109.471
angle_coeff 11 2.168 109.471
angle_coeff 12 2.168 109.471
angle_coeff 13 2.168 109.471
angle_coeff 14 2.168 106.700
angle_coeff 15 2.168 106.700
angle_coeff 16 2.168 104.510
angle_coeff 17 2.168 109.471

dihedral_style harmonic
dihedral_coeff 1 0.9756 -1 2
dihedral_coeff 2 0.9756 -1 2
dihedral_coeff 3 0.9756 -1 2
dihedral_coeff 4 0.9756 -1 2
dihedral_coeff 5 0.9756 -1 2
dihedral_coeff 6 0.02168 -1 6
dihedral_coeff 7 0.02168 -1 6
dihedral_coeff 8 0.00000 -1 6
dihedral_coeff 9 0.00000 -1 3
dihedral_coeff 10 0.0433 -1 3
dihedral_coeff 11 0.0433 -1 3
dihedral_coeff 12 0.9756 -1 2
dihedral_coeff 13 0.9756 -1 2
dihedral_coeff 14 0.9756 -1 2
dihedral_coeff 15 0.9756 -1 2
dihedral_coeff 16 0.02168 -1 6
dihedral_coeff 17 0.02168 -1 6
dihedral_coeff 18 0.00000 -1 6
dihedral_coeff 19 0.00000 -1 6
dihedral_coeff 20 0.04336 -1 3
dihedral_coeff 21 0.9756 -1 2
dihedral_coeff 22 0.9756 -1 2
dihedral_coeff 23 0.9756 -1 2
dihedral_coeff 24 0.9756 -1 2
dihedral_coeff 25 0.02168 -1 6
dihedral_coeff 26 0.02168 -1 6

velocity fluid create 300.0 1929457 rot yes dist gaussian

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes page 500000 one 50000

fix pistonfreeze piston setforce 0.0 0.0 0.0
fix membranefreeze membrane spring/self 1000
compute selective_thermostat thermostat_target temp
thermo_modify temp selective_thermostat

thermo_style one
thermo 10
run 0

minimize 1.0e-4 1.0e-6 1000 1000

fix relax all box/relax x 0.0 y 0.0
minimize 1.0e-4 1.0e-6 100 1000

unfix relax
unfix pistonfreeze

fix pistonkeep piston setforce 0.0 0.0 NULL
fix piston1thrust piston1 aveforce NULL NULL 0.0005
fix piston2thrust piston2 aveforce NULL NULL -0.0005
fix 1 fluid shake 1.0e-4 100 0 b 7 a 16
fix NVTequilib thermostat_target nvt temp 300 300 0.05

timestep 0.0005
thermo 1000
dump 1 all xyz 1000 dump1.xyz
run 100000

unfix NVTequilib
unfix piston1thrust
unfix piston2thrust

reset_timestep 0
timestep 0.0005
fix piston1push piston1 aveforce NULL NULL 0.000769
fix piston2push piston2 aveforce NULL NULL -0.0000015
fix NVT thermostat_target nvt temp 300 300 0.05

region W block -200 200 -200 200 1 200 units box
region Q block -200 200 -200 200 -200 -1 units box
group Water1 dynamic H2O region W
group Water2 dynamic H2O region Q
group Ion dynamic salt region W
variable Water1Number equal count(Water1)
variable IonNumber equal count(Ion)
variable Water2Number equal count(Water2)

akohlmey · February 28, 2022, 9:59pm

so and what do you expect that should happen now?

majido · March 1, 2022, 6:23am

Dear Axel,

I expected the water molecules to be pressed against the membrane surface by the piston (500 bar pressure). But they do not spread in the simulation box and resist the piston like a column as I attach Fig.

tnk

akohlmey · March 1, 2022, 6:47am

That is not really answering my question. Let me rephrase:

What kind of response do you expect from the mailing list?

To summarize:

In your original post you just state that your water molecules “shrink” despite using fix shake, and then claim they are “resisting” a “piston pressure”, but do not ask any questions and instead quote a lengthy and confusing input file without any explanations or any way to test it
Now you are producing a picture without any explanation of what is represented by it and how it corresponds to your previous post. And again there are no questions or hints at what you expect to get from the mailing list.

So how can anybody make sense of that? or figure out what you want?
Obviously you are not happy with your simulation, but do you expect people now not only to provide answers, but also to supply the questions?