[ESPResSo-users] NVT with iccp3m

[Xikai Jiang]

Hello:

Sorry to bother who are not interested. I'm trying to do a NVT simulation of VdW liquid between two dielectric plates using iccp3m. While running the simulation, the output total energy is NaN and the temperature is higher than what I want (I set T=1 for every atom between two walls, LANGEVIN_PER_PARTICLE is enabled). Below is a portion of the output:

t=6.399999999998423 E=-NaN, T=71.43363446128329

t=6.49999999999819 E=-NaN, T=69.94231211553624

t=6.599999999997957 E=-NaN, T=60.330816719961994

Previously when I set T=450, the simulation couldn't run with a error of segmentation fault. Then I set T=1, the simulation can run, but the temperature is not controlled well. Not sure where I went wrong.

I'm using the latest development code that is downloaded from Github. I appreciate any inputs. Thanks!  - Jimmy

Please find the script for my simulation below:

# units:

# length (nm)

# time (ps)

# energy (KJ/mol)

source "tests_common.tcl"

# basic settings

set box_lx 3.; set box_ly 3.; set box_lz 5.

set spacing 0.15; set Natom  [expr int($box_lx/$spacing)+1]

# number of total wall atoms:

set Nstart [expr $Natom*$Natom*2]

set density 0.3

# number of all charged particles:

set total_ions [expr 2*ceil($box_lx*$box_ly*2.9*$density)]

set epsilonr 2

setmd box_l $box_lx $box_ly $box_lz

set dt 0.0001; setmd time_step $dt

setmd skin 0.3

set temp 1; set gamma 1

thermostat langevin $temp $gamma

if { [regexp "ROTATION" [code_info]] } {

    set deg_free 6

} {

    set deg_free 3

}

# set up walls

set counter 0

set normals [ list ]; set areas [ list ]

set epsilons [ list ]; set ext_fields [ list ]

#bottom wall

for { set i 0 } { $i < $Natom } { incr i } {

    for { set j 0 } { $j < $Natom } { incr j } {

part $counter pos [expr $i*$spacing] [expr $j*$spacing] 0 fix 1 1 1 q 0.1 type 0

lappend normals { 0 0 1. }

        lappend areas [expr $spacing*$spacing]

        lappend epsilons 10000

        incr counter

    }

}

#top wall

for { set i 0 } { $i < $Natom } { incr i } {

    for { set j 0 } { $j < $Natom } { incr j } {

        part $counter pos [expr $i*$spacing] [expr $j*$spacing] 2.9 fix 1 1 1 q 0.1 type 0

        lappend normals { 0 0 -1. }

        lappend areas [expr $spacing*$spacing]

        lappend epsilons 10000

        incr counter

    }

}

# purely repulsive charges

set q 1; set type 1

for {set i 0} {$i < $total_ions} {incr i} {

    set posx [expr $box_lx*[t_random]]

    set posy [expr $box_ly*[t_random]]

    set posz [expr (2.9-1)*[t_random]+0.5]

    set q [expr -$q];set type [expr 1-$type]

    # temperature for each charged atom between two walls

    part [expr $Nstart+$i] pos $posx $posy $posz  q $q temp $temp gamma $gamma type [expr $type+1]

}

# LJ interactions

set sig 1.; set cut [expr 1.12246*$sig]

set eps 1

inter 0 1 lennard-jones $eps $sig $cut

inter 0 2 lennard-jones $eps $sig $cut

inter 1 1 lennard-jones $eps $sig $cut

inter 1 2 lennard-jones $eps $sig $cut

inter 2 2 lennard-jones $eps $sig $cut

# forcecapping

set integ_steps 1000

inter forcecap individual

for {set i 1} {$i < 10} {incr i} {

    puts "entered forcecap loop"

    set rad [expr 1.0 - 0.5*$i/10.0]

    set lb [expr 138.935/$epsilonr*$i/10.0]

    inter 1 1 lennard-jones $eps $sig $cut 0 0 $rad

    inter 1 2 lennard-jones $eps $sig $cut 0 0 $rad

    inter 2 2 lennard-jones $eps $sig $cut 0 0 $rad

    puts "[inter coulomb [expr $lb/$temp] p3m tunev2 accuracy 1e-3 r_cut 0 mesh 0 cao 0]"

    integrate $integ_steps

    puts "forcecap $i,time:[expr $i*$integ_steps*$dt]ps"

}

inter forcecap 0

puts "[inter coulomb [expr 138.935/$temp/$epsilonr] p3m tunev2 accuracy 1e-3 r_cut 0 mesh 0 cao 0]"

# ICCP3M, zero potential on each wall

iccp3m $Nstart eps_out 1 max_iterations 60 convergence 1e-3 relax 0.7 areas $areas normals $normals epsilons $epsilons

# MD integration

for {set i 0} { $i < 100 } { incr i} {

    set temp [expr [analyze energy kinetic]/(($deg_free/2.0)*$total_ions)]

    puts "t=[setmd time] E=[analyze energy total], T=$temp"

    integrate $integ_steps

}

exit 0