|
From: | Xikai Jiang |
Subject: | [ESPResSo-users] NVT with iccp3m |
Date: | Mon, 17 Mar 2014 15:19:52 -0400 |
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 |
[Prev in Thread] | Current Thread | [Next in Thread] |