Re: problems with simulation of LJ fluid

From: Axel Kohlmeyer (akohlmey_at_gmail.com)
Date: Mon Dec 03 2012 - 11:20:58 CST

On Sun, Dec 2, 2012 at 12:30 AM, <wenglindong_at_yahoo.com.cn> wrote:
> Hi, NAMDers,
>
> Is there anyone who ever simulated Lennard Jones fluids, especially pure Ar
> or pure Kr, using NAMD? I tried to do so, but got wrong results for Ar/Kr
> mixtures and even could not run the equilibrium of pure Kr or pure Ar
> fluids. Attached is the top, par, conf files, any suggestions would be
> really appreciated.

your starting geometry is bad.
you are likely to have some very close contacts.
most likely due to incorrect choice of cell dimensions.

you can't blame NAMD for that, if you feed it bad
data for input, i can't make it right.

axel.

>
> Topology file:
>
> !Ar/Kr
> MASS 1 AR 39.948000 AR !
> MASS 2 KR 83.798000 KR !
>
> DEFA FIRS none LAST none
> AUTOGENERATE ANGLES DIHEDRALS
>
> RESI ARLJ 0.00 ! Ar Atom
> GROUP
> ATOM AR1 AR 0.00
> PATCHING FIRST NONE LAST NONE
>
> RESI KRLJ 0.00 ! Kr Atom
> GROUP
> ATOM KR1 KR 0.00
> PATCHING FIRST NONE LAST NONE
>
> END
>
> Parameter file:
> NONBONDED nbxmod 5 atom cdiel shift vatom vdistance vswitch -
> cutnb 14.0 ctofnb 12.0 ctonnb 10.0 eps 1.0 e14fac 1.0 wmin 1.5
> !adm jr., 5/08/91, suggested cutoff scheme
> !
> !V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]
> !
> !epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j)
> !Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j
> !
> !atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
> !
> ! ions, note lack of NBFIXes
> !
> AR 0.0 -0.2381 1.91099 ! Ar
> KR 0.0 -0.3319 2.03895 ! Kr
>
>
> HBOND CUTHB 0.5 ! If you want to do hbond analysis (only), then use
> ! READ PARAM APPEND CARD
> ! to append hbond parameters from the file: par_hbond.inp
>
> END
>
> Config file:
> #############################################################
> ## JOB DESCRIPTION ##
> #############################################################
>
> # Minimization and Equilibration of a Kr Box
>
> #############################################################
> ## ADJUSTABLE PARAMETERS ##
> #############################################################
>
> structure kr_b.psf
> coordinates kr_b.pdb
>
> set temperature 94.4 ;# K
> set outputname kr_b
>
> #############################################################
> ## SIMULATION PARAMETERS ##
> #############################################################
>
> temperature $temperature ;# initialize velocities randomly
> firsttimestep 0
>
> # Force-Field Parameters
> paraTypeCharmm on
> parameters par_lj.prm
> exclude scaled1-4
> 1-4scaling 1.0
>
> # You have some freedom choosing the cutoff
> cutoff 12. ;# may use smaller, maybe 10, with PME
> switching on
> switchdist 10. ;# cutoff - 2
> pairlistdist 13.5 ;# > cutoff
>
> # Integrator Parameters
> timestep 2.0 ;# 2 fs/step
> rigidBonds all ;# needed for 2 fs steps
> nonbondedFreq 1 ;# nonbonded forces every step (>= 2 fs, here one
> step)
> fullElectFrequency 2 ;# PME only every other step
> stepspercycle 10 ;# redo pairlists every ten steps
>
> # Constant Temperature Control
> langevin on ;# langevin dynamics
> langevinTemp $temperature ;# random noise at this level
> langevinDamping 5. ;# damping coefficient of 5 1/ps
> langevinHydrogen off ;# don't couple langevin bath to hydrogens
>
> # Constant Pressure Control
> useGroupPressure yes ;# require periodic boundary conditions
> useFlexibleCell no ;# no for isotropic system, yes for anisotropic one
> useConstantRatio no ;# no for isotropic system, yes for anisotropic one
> useConstantArea no ;# no for isotropic system, yes for anisotropic one
>
> # (Should be accompanied by Langevin dynamics for temperature control)
> LangevinPiston on
> LangevinPistonTarget 1.01325
> LangevinPistonPeriod 100. ;# oscillation period around 100 fs
> LangevinPistonDecay 50. ;# oscillation decay time of 50 fs
> LangevinPistonTemp $temperature ;# equal to the controlled temperature
>
> # Periodic Boundary Conditions
> cellBasisVector1 32. 0. 0.
> cellBasisVector2 0. 32. 0.
> cellBasisVector3 0. 0. 32.
> cellOrigin 0. 0. 0.
> wrapAll on ;# wrap all molecules
>
> # PME (for full-system periodic electrostatics)
> PME yes
> PMEGridSizeX 32
> PMEGridSizeY 32
> PMEGridSizeZ 32
>
> # Output
> outputName $outputname
> restartfreq 1000 ;# 1000 steps = every 2 ps
> dcdfreq 10
> velDCDfreq 10
> xstFreq 500
> outputEnergies 500 ;# 500 steps = every 1 ps
> outputTiming 500 ;# shows time per step and time to completion
> outputPressure 500
>
> #############################################################
> ## EXECUTION SCRIPT ##
> #############################################################
>
> # Minimization
> minimize 10000 ;# lower potential energy for 10000
> steps (20 ps)
> reinitvels $temperature ;# since minimization zeros velocities
>
> run 500000 ;# 1 ns
>
> The wrong information about the simulation of pure Kr fluids is:
> ... ...
> LINE MINIMIZER BRACKET: DX 4.08248e-302 3.67423e-301 DU 0 0 DUDX
> -2.21951e+11 -2.21951e+11 -2.21951e+11
> LINE MINIMIZER REDUCING GRADIENT FROM 3.32927e+11 TO 3.32927e+08
> PRESSURE: 10000 1.00128e+06 7032.42 -16055.9 7032.42 1.37011e+06 -28789.2
> -16055.9 -28789.2 4.01709e+06
> GPRESSURE: 10000 1.00128e+06 7032.42 -16055.9 7032.42 1.37011e+06 -28789.2
> -16055.9 -28789.2 4.01709e+06
> TIMING: 10000 CPU: 655.735, 0.0628918/step Wall: 655.735, 0.0628918/step,
> 0 hours remaining, 726.261719 MB of memory in use.
> ETITLE: TS BOND ANGLE DIHED IMPRP
> ELECT VDW BOUNDARY MISC KINETIC
> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
> PRESSURE GPRESSURE VOLUME PRESSAVG GPRESSAVG
>
> ENERGY: 10000 0.0000 0.0000 0.0000 0.0000
> 0.0000 9999999999.9999 0.0000 0.0000 0.0000
> 9999999999.9999 0.0000 9999999999.9999 9999999999.9999
> 0.0000 2129495.3815 2129495.3815 32768.0000 2129495.3815
> 2129495.3815
>
> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 10000
> WRITING COORDINATES TO DCD FILE AT STEP 10000
> WRITING COORDINATES TO RESTART FILE AT STEP 10000
> FINISHED WRITING RESTART COORDINATES
> WRITING VELOCITIES TO DCD FILE AT STEP 10000
> WRITING VELOCITIES TO RESTART FILE AT STEP 10000
> FINISHED WRITING RESTART VELOCITIES
> REINITIALIZING VELOCITIES AT STEP 10000 TO 94.4 KELVIN.
> TCL: Running for 500000 steps
> PRESSURE: 10000 1.03964e+06 6761.74 -21242.5 6761.74 1.46064e+06 -106189
> -21242.5 -106189 5.06431e+06
> GPRESSURE: 10000 1.03964e+06 6761.74 -21242.5 6761.74 1.46064e+06 -106189
> -21242.5 -106189 5.06431e+06
> ETITLE: TS BOND ANGLE DIHED IMPRP
> ELECT VDW BOUNDARY MISC KINETIC
> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
> PRESSURE GPRESSURE VOLUME PRESSAVG GPRESSAVG
>
> ENERGY: 10000 0.0000 0.0000 0.0000 0.0000
> 0.0000 9999999999.9999 0.0000 0.0000 206.8176
> 9999999999.9999 95.1763 9999999999.9999 9999999999.9999
> 95.1763 2521528.8179 2521528.8179 32768.0000 2521528.8179
> 2521528.8179
>
> ERROR: Atom 12 velocity is 1.03923e+12 -3.01387e+11 -7.6917e+11 (limit is
> 6000, atom 0 of 257 on patch 0 pe 8)
> ERROR: Atom 28 velocity is 126.522 -1146.48 12682.8 (limit is 6000, atom 3
> of 257 on patch 0 pe 8)
> ERROR: Atom 55 velocity is -1.50636e+11 -2.28264e+11 3.73979e+11 (limit is
> 6000, atom 10 of 257 on patch 0 pe 8)
> ERROR: Atom 180 velocity is 10799 483020 1.23976e+06 (limit is 6000, atom 35
> of 257 on patch 0 pe 8)
> ... ...
> ERROR: Atom 19 velocity is -1.27025e+11 6.85185e+11 6.58609e+11 (limit is
> 6000, atom 223 of 232 on patch 2 pe 4)
> ERROR: Atom 27 velocity is 1.27025e+11 -6.85185e+11 -6.58609e+11 (limit is
> 6000, atom 224 of 232 on patch 2 pe 4)
> ERROR: Atom 109 velocity is -1.01188e+12 1.09255e+12 -2.67297e+11 (limit is
> 6000, atom 226 of 232 on patch 2 pe 4)
> ERROR: Atom 117 velocity is 1.01188e+12 -1.09255e+12 2.67297e+11 (limit is
> 6000, atom 227 of 232 on patch 2 pe 4)
> ERROR: Atom 684 velocity is -110.949 -691.142 -13593.1 (limit is 6000, atom
> 229 of 232 on patch 2 pe 4)
> ERROR: Atoms moving too fast; simulation has become unstable (62 atoms on
> patch 2 pe 4).
> ERROR: Exiting prematurely; see error messages above.
> ====================================================
>
>
> Lindong
> COE
> University of North Carolina at Charlotte
>
>

--
Dr. Axel Kohlmeyer  akohlmey_at_gmail.com  http://goo.gl/1wk0
International Centre for Theoretical Physics, Trieste. Italy.

This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:22:49 CST