LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Density of methanol

# Re: [lammps-users] Density of methanol

 From: Axel Kohlmeyer Date: Thu, 17 Aug 2017 10:45:58 -0400

On Thu, Aug 17, 2017 at 10:15 AM, 노승효 wrote:

To lammps users,

I am testing the density of methanol using a given example of LAMMPS software (~/LAMMPS/Examples/dreading/)

It is known that the density of methanol is ~0.792 g/cm3.

The initial density of methanol example was 0.457 g/cm3 but the density is lower, showing the enlarged unitcell.

How can I get reasonable vlaue (~0.792 g/cm3)?

​different parameterizations in classical MD will reproduce experimentally data differently well.
there is no compound and not set of force field parameters that reproduces every experimental property exactly. parameterizations are always a compromise.
in simple empirical potentials, the parameterization is often limited as to what themodynamic states can be reproduced well.

to find a parameterization that reproduces methanol more to your liking, you will have to study the published literature. folks have been studying methanol (and methanol mixtures) for as long as i am doing simulations, i.e. over 20 years. so there should be plenty of literature to check out.

​axel.​

## results##

Step Temp TotEng PotEng KinEng Lx Ly Lz Press Density
0       298.15    12308.081    3094.6323    9213.4485     59.99067     57.38448     58.40913    1781.0418    0.4572555
1000    290.58157    19927.153    10947.585    8979.5683    70.096033    67.050833    68.248084    -92.74926   0.28663535
2000    305.64025    21030.387    11585.475    9444.9127    81.307743     77.77547    79.164219    60.844666    0.1836599
3000    302.94026     21826.71    12465.233    9361.4774    93.059138    89.016346    90.605811    97.386573   0.12249913

## input script ##

units           real
atom_style      full
boundary        p p p
dielectric      1
special_bonds   lj/coul 0.0 0.0 1.0

pair_style      hybrid/overlay hbond/dreiding/lj 2 6 6.5 90 lj/cut/coul/long 8.50000  11.5
bond_style      harmonic
angle_style     harmonic
dihedral_style  harmonic
improper_style  none
kspace_style    pppm 0.001

replicate 3 3 3

pair_coeff      1    1    lj/cut/coul/long        0.015200000256300         2.846421344984478
pair_coeff      1    2    lj/cut/coul/long        0.001232882795416         2.846421344984478
pair_coeff      1    3    lj/cut/coul/long        0.038019995160237         3.159705878878677
pair_coeff      1    4    lj/cut/coul/long        0.038139744011598         2.939787518071103
pair_coeff      2    2    lj/cut/coul/long     9.99999974737875e-05         2.846421344984478
pair_coeff      2    3    lj/cut/coul/long        0.003083828758188         3.159705878878677
pair_coeff      2    4    lj/cut/coul/long        0.003093541672406         2.939787518071103
pair_coeff      3    3    lj/cut/coul/long        0.095100000500679         3.472990412772877
pair_coeff      3    4    lj/cut/coul/long        0.095399530150179         3.253072051965302
pair_coeff      4    4    lj/cut/coul/long        0.095700003206730         3.033153691157727
pair_coeff      4    4    hbond/dreiding/lj 2 i 0.4000E+01 2.750000000000000 4
pair_modify     mix arithmetic
neighbor        2.0 multi
neigh_modify    every 2 delay 4 check yes
variable        input index in.ch3oh.box.dreiding
variable        sname index ch3oh.box.dreiding

compute   hb all pair hbond/dreiding/lj
variable    C_hbond equal c_hb[1] #number hbonds
variable    E_hbond equal c_hb[2] #hbond energy
#thermo_style   custom etotal ke temp pe ebond eangle edihed eimp evdwl ecoul elong v_E_hbond v_C_hbond press vol
#thermo_modify  line multi format float %14.6f

# Output
thermo         1000
thermo_style   custom step temp etotal pe ke lx ly lz press density

dump myOvito all custom 1000 dump.ovito.* id type xs ys zs

################### NPT
timestep 1
velocity       all create 298.15 4928459 rot yes mom yes dist gaussian
#fix ensemble_press all press/berendsen iso 1.0 1.0 1000.0
fix   fxnpt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 100.0
run 100000
unfix fxnpt
#unfix ensemble_press
reset_timestep 0

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most