Trouble to find water diffusion rate

Continuing the discussion from Trouble finding water diffusion:

#Input file for chloride diffusion in jennite crystal using MSD method

units real
atom_style atomic
boundary p p p
read_data jennitecombined3.lmp
neigh_modify delay 0 every 1

Define the group for chloride atoms

group chloride type 8
group csh1 id 1 2 3 4 1:690
group water id 5 6 691:1770
group csh2 id 1 2 3 4 1771:2460
group sodiumchloride id 7 8 2461:2466

region 1 block 0.0 28.00 -2.00 35.00 -3.0 10.85
group csh1 region 1
region 2 block 0.0 28.00 -2.00 35.00 10.85 49.80
group water region 2
region 3 block 0.0 28.00 -2.00 35.00 49.80 60.00
group csh2 region 3

Fix the positions of groups csh1 and csh2

fix fix_csh1 csh1 setforce 0.0 0.0 0.0
fix fix_csh2 csh2 setforce 0.0 0.0 0.0
#fix fix_water water setforce 0.3 0.0 0.0

Set the communication cutoff to 10 Angstrom

comm_modify cutoff 5

#Define the force field
pair_style lj/cut 2.5

neighbor 2.5 bin

pair_coeff * * 0.294397063 0.033562214
pair_coeff * * 0.38501583 0.033692483
pair_coeff * * 0.287659646 0.041740485
pair_coeff * * 0.407152299 0.002963415
pair_coeff * * 0.497972184 0.002967482
pair_coeff * * 0.378270062 0.521398045
pair_coeff * * 0.316628818 0.650037879
pair_coeff * * 0.436240376 0.003699858
pair_coeff * * 0.0 0.0
pair_coeff * * 0.38501583 0.0
pair_coeff * * 0.436240376 0.003699858
pair_coeff * * 0.316628818 0.650037879
pair_coeff * * 0.323239045 0.042082192
pair_coeff * * 0.323239045 0.042082192
pair_coeff * * 0.294397063 0.0
pair_coeff * * 0.442973487 0.000238235
pair_coeff * * 0.330362359 0.002692308
pair_coeff * * 0.55591693 2.10055E-05
pair_coeff * * 0.378270062 0.521398045
pair_coeff * * 0.31642861 0.652509506
pair_coeff * * 0.294397063 0.0
pair_coeff * * 0.38501583 0.0
pair_coeff * * 0.258328345 0.418925676
pair_coeff * * 0.349474969 0.41580163
pair_coeff * * 0.287402102 0.523458647
pair_coeff * * 0.38501583 0.0
pair_coeff * * 0.294397063 0.0
pair_coeff * * 0.44041546 0.414524349

set time step

timestep 0.005 # time step of 1 fs

#Set up nvt ensemble
fix 1 all nvt temp 300 300 100

velocity chloride create 300 1234 mom yes dist uniform

Set up the minimization settings

min_style cg
min_modify line quadratic

Run energy minimization

minimize 1.0e-10 1.0e-10 1000 10000

Equilibrate the system for 50000 time steps

thermo 1000
run 500000

#Data gathering steps
reset_timestep 0

Set up the compute commands

compute msd all msd com yes

Set up the fix and compute commands

compute msd_water water msd com yes
variable twopoint equal c_msd[4]/4/(stepdt+1.0e-6)
fix 9 all vector 10 c_msd[4]
variable fitslope equal slope(f_9)/4/(10

thermo_style custom step temp c_msd[4] v_twopoint v_fitslope

Set up the dump file for all atoms

dump 1 all atom 100

#Set up the dump file for chloride atoms
dump 2 water atom 100

Run the simulation

thermo 1000
run 1000000

This code is calculating all the diffusion in same amount. I mean msd of chloride ,water and all other particles are showing same result.

That’s a weird script. To pick just a few points:

  • You are applying each pair_coeff to all atom types by using the starts * *, which I never seen before.
  • You only have one single compute msd applied to the group water, but your message refers to the msd of all the other species as well.
  • The group definitions are suspicious… are you water molecule ids really 5 6 691:1770 ?

My advice, take a step back, follow tutorials and try to understand how the basics of LAMMPS work.

1 Like

@Noshin Please do not open a new thread related to the same issue.

You already opened a discussion, consider sticking to your initial thread to avoid useless multiplication of forum topics.

These are pointless. Each command with pair_coeff * * will wipe out the settings from the previous one.

Why do you need this?

What kind of water model only needs atom style “atomic”, i.e. has no bonds?

You say you are also modeling ions, but use a model without charges. This is highly unusual.

This is a very short cutoff for real units.

These are very unusual values for epsilon and sigma for an atomic/molecular model in real units.

1 Like

Dear simongravelle,
In my data file the water molecules are defined as H and O. I gave them different id so they dont overlap with the H and O of my crystal and yes, the group water is in 691:1770.
I am new to this field and I was doing this for my undergrad project, I decided to choose the Atomic style atomic and force field lj to keep it simple at first. I wrote this script by reading the manual book of lammps . As I am new, I may mistake to understand the pair_coeff part and thank you for mentioning me.

Thank you for your valuable suggestions. I am currently doing the Udemy course for lammps.

Dear akohlmey,
You have pointed out my errors and thank you for that.
I did not understand the pair_coeff at first. I am studying this and correcting this point.
At first I wrote my script without comm_modify cutoff and my code didnot run, on the ouput screen there was mentioned that I must have introduced comm_modify.
I’m new to this field. So, I wanted to run a basic and simple model at first. As I am new, I may mistake and if I start from the basic I can solve it fast.
Thank you for your valuable suggestions. I am trying to correct my existing script.

Dear Germain,
I am extremely sorry for my inconvenience. I tried to undo my previously opened thread but unfortunately it’s not deleting. I will make sure not to make further trouble.

You have picked a very difficult project for an undergraduate with no prior experience. This is not your fault – your course instructors should not let you wander into the wild forest alone (after all, this is what you pay them course fees for). The biggest problem you will face is that, since diffusion inside a solid is very slow, it is impossible for you to get sensible results without access to supercomputing resources.

I would suggest an alternative project. This paper: studies molecular dynamics simulations of ionic conductivity in salt water. The authors have made their LAMMPS scripts publicly available. Unfortunately those scripts use custom modified commands for the analysis – but the simulations themselves are sound. You could write a simple but useful undergraduate project along the following lines:

  1. The initial configurations are for modelling salt water at concentrations of 0.5, 1.0, 2.5 and 4.0 M. Look through the data files and, counting the number of ions and dividing by the box volume, confirm that the concentrations are accurate.

  2. Starting from the supplied scripts, write a script that runs a simulation for 10 nanoseconds at a temperature of 300 K and a pressure of 1 bar. Underneath each line of your script, write a comment explaining what it does, including the effect of changing any numbers.

  3. Run a set of experiments that changes one “variable” and studies the effect on density – for example density vs temperature at fixed pressure and fixed salt concentration, or density vs salt concentration at fixed pressure and temperature. Compare your results to experimental values of the density of salt water. Do you get a similar trend?

Dear stree,
Thank you soo much. Your advice is very helpful. I’m really glad you explained these points so well and in a easy way for me. I’m studying on this new topic.