Re: problems with simulation of LJ fluid

From: Rajan Vatassery (rajan_at_umn.edu)
Date: Mon Dec 03 2012 - 12:35:42 CST

Lindong,
        Minimizing 10000 steps is a really long time and almost certainly
should have reduced the close contacts effectively. You probably have
some superimposed atoms (atoms that are exactly the same coordinates)
because close contacts wouldn't make the large energies you're finding.
A very poor (small) choice of periodic box might possibly be the
culprit, but try to look through your initial input for superimposed
atoms. VMD should be very useful in this respect.

Good luck,

rajan

On Sun, 2012-12-02 at 07:30 +0800, 翁林岽 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.
>
>
> 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
>
>
>
>

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