LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Water diffusion coefficient from MSD
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Water diffusion coefficient from MSD


From: Ray Shan <rshan@...1795...>
Date: Thu, 19 Oct 2017 01:12:38 +0000

Perhaps you forgot to add fix nve to your script.

 

Ray

 

From: "sungpar@...4084..." <sungpar@...4084...>
Date: Wednesday, October 18, 2017 at 3:06 AM
To: Ray Shan <rshan@...1795...>
Subject: RE: [lammps-users] Water diffusion coefficient from MSD

 

Dear,

Thanks for your advises.

First, I ran the input file with npt 1ns to equilibriate the system, final density was 0.99 g/cc. So, now I think that system is stable.

And then, I ran the file with 1, nve with temp/restart 2, nvt with deform, 3 npt for 300ps each.

MSD results are almost same, ther result within 200ps MSD should be around 200-300 (angstrom^2).

But, my result is about 1-2 (angstrom^2).

Specially, when I ran with nve, the calculated kinetic energy is not changed during calculation, and I think it means the molecular doesnt move I thing I did something very wrong. Hoding this problem almost 2 weeks, can you give me any more advises ?

 

Thanks,

 

 

units                           real

boundary                     p p p

 

atom_style                    full              

pair_style                      lj/cut/coul/cut 10 12

bond_style              harmonic

angle_style              harmonic

 

read_restart                  restart100W.txt

 

neighbor                      2.0 bin

neigh_modify          every 10 delay 0 check yes

 

thermo_style                 custom step pe ke etotal temp press vol density

thermo           100

 

timestep          1

 

fix                         1 all shake 1.0e-4 10 0 b 1 a 1

 

#fix                         2 all deform 1 x final 0.0 14.937 y final 0.0 14.937 z final 0.0 14.937

fix                               3 all npt temp 298.15 298.15 100 iso 1 1 100

#fix                             3 all nvt temp 298.15 298.15 100

#fix                          3 all nve

#fix                          10 all temp/rescale 100 298.15 298.15 0.05 1.0

 

#run     1000

#unfix   2

#run     1000

 

compute                     1 all msd

 

fix            16 all ave/time 1 100  1000 c_1[1] c_1[2] c_1[3] c_1[4] file msd.txt

 

dump            1 all xyz 100 dump.xyz

dump_modify 1 element O H C

 

run                  300000

 

 

 

 

 

 

보낸 사람: Ray Shan
보낸 날짜: 2017 10 17일 화요일 오전 11:44
받는 사람: sungpar@...4084...
참조: lammps-users@lists.sourceforge.net
제목: Re: [lammps-users] Water diffusion coefficient from MSD

 

Many things have gone wrong, to name a few:

 

1.     A density of 0.9 g/cc indicates the inadequacy of the chosen potential for describing liquid water molecules.

2.     100 ps of NPT dynamics is too short.

3.     MSD for diffusion coefficient should be ran with NVE ensemble with the thermostat.

4.     Diffusion run of 100-200 ps is far too short.

 

Ray

 

From: "sungpar@...4084..." <sungpar@...4084...>
Date: Monday, October 16, 2017 at 4:32 AM
To: "lammps-users@lists.sourceforge.net" <lammps-users@lists.sourceforge.net>
Subject: [lammps-users] Water diffusion coefficient from MSD

 

Dear, Lammps uers

 

I am trying to calculate self diffusion coefficient from MSD.

 

I made 100 water molecular(liquid) 1 atm, 25 C.

First of all, I gave time 100ps to equilibriate the system with npt.

I got density about 0.9 g/ml3, and 3332 angstrom^2 volume.

 

And next, I simulate the system using nvt. I fixed the system volume with fix deform for 10000 run.

 

Then, I used msd command with O atoms.

But, diffusion coefficient is too low, it is 1/100 than I expected.

I ran the simulation 200ps. ( I tried much more time, but result was same)

 

Can any one give advises ?

 

### Units

units                           real

### Resion

boundary                     p p p

### Force type

atom_style                    full              

pair_style                      lj/cut/coul/cut 9 9

bond_style              harmonic

angle_style              harmonic

### Atom definition

read_data                     data_100W.txt

velocity             all create 300 8343 dist gaussian units box

### Neighbor list

neighbor                      2.0 bin

neigh_modify          every 10 delay 0 check yes

### Calculation

thermo_style                 custom step pe ke etotal temp press vol density

thermo           10

### Minimization

minimize                      1.0e-5 1.0e-7 1000 100000

### Timestep

timestep          1

### Fix molecular

fix                         1 all shake 1.0e-4 10 0 b 1 a 1

### Dynamic type

fix                               2 all npt temp 298.15 298.15 100 iso 1 1 100

### Variable

variable fs equal 1*step

variable density equal "density"

variable vol equal "vol"

### Print density

fix                         3 all print 10 "${fs} ${density}"    file fsDensity.txt         screen no

fix                         4 all print 10 "${fs} ${vol}"    file fsVol.txt         screen no

### Start running

run                  100000

### Save the restart file

write_restart restart100W.txt

 

 

### Units

units                           real

### Resion

boundary                     p p p

### Force type

atom_style                    full              

pair_style                      lj/cut/coul/cut 9 12

bond_style              harmonic

angle_style              harmonic

### Atom definition

read_restart                  restart100W.txt

### Neighbor list

neighbor                      2.0 bin

neigh_modify          every 10 delay 0 check yes

### Calculation

thermo_style                 custom step pe ke etotal temp press vol density lx ly lz

thermo           100

### Timestep

timestep          1

### Fix molecular

fix                         1 all shake 1.0e-4 10 0 b 1 a 1

### Dynamic type

 

fix                         2 all deform 1 x final 0.0 14.937 y final 0.0 14.937 z final 0.0 14.937

fix                               3 all nvt temp 298.15 298.15 100

 

run     10000

unfix     2

run     10000

 

### Group

group                    ox type 1

### Compute

compute                     1 ox msd

compute                 5 ox vacf

### Variable

variable  msd            equal c_1[4]

variable  vacf            equal c_5[4]

variable  fs               equal 1*step

variable  Temp          equal "temp"

variable  Vol              equal "vol"

variable  density         equal "density"

### Print results and make files

fix          7  all print 10  "${fs} ${msd}"         file msd.txt         screen no

fix            9  all print 10  "${fs} ${vacf}"     file vacf.txt         screen no

fix            8  all vector 1 c_5[4]

variable    diff equal dt*trap(f_8)

fix            10 all print 10 "${fs} ${diff}"     file diff.txt         screen no

 

### Dump

dump            1 all xyz 100 dump.xyz

dump_modify 1 element O H C

### Run a simulation

run                  200000

 

 

id:image001.png@...7171...