[ESPResSo-users] crash using RATTLE with MPI

[Muhammad Anwar]

Hello, I am running espresso-3.0.1 with lam/MPI.

All tests pass both in serial and parallel, but when I try to run a 
calculation with restrained bonds on more than 1 cpu, it breaks with the 
output at the end of this email.

The important lines of the script are:

set bond_l [expr 1.53]
inter 1 rigid_bond $bond_l 1e-2 1e-2

If I run on 1 CPU, output is OK:
.
.
1.530076435924329 Warmup finished. Minimal distance now 1.530076435924329
Step 609/200^CCtrl-C caught...


If I run with the rigid_bond replaced with a harmonic one, output is OK 
in parallel or serial:
inter 1 harmonic 1 7
.
.
.
2.738430991400251 Warmup finished. Minimal distance now 2.738430991400251
Step 6061/20000000      6.06000 -0.0082 66.4397  6602.3605  65.5916    
-0.253^CCtrl-C caught...

Have I found a bug, or am I doing something wrong?
Thanks
Best Regards,
Muhammad Anwar

----------------------------------------

Interactions:
  {1 RIGID_BOND 1.53 0.01 0.01} {2 angle 134.19148789844388 
1.9163715186897738} {3 dihedral 1 -0.969533500066257 3.623170173257985} 
{0 0 lj-gen 0.10449043857692163 4.5 9.0 0.031 0.0 12  6  1.0 2.0 0.0 } 
{0 1 lj-gen 0.1627071790768632 4.5 9.0 0.031 0.0 12  6  1.0 2.0 0.0 } {1 
1 lj-gen 0.25317460716839746 4.5 9.0 0.031 0.0 12  6  1.0 2.0 0.0 }
Creating vmd_start.script...
vmd_start.script finished.
background_errors 0 {053 RATTLE failed to converge after 1000 
iterations} {004 some particles moved more than min_local_box_l, reduce 
the time step} 1 <consent> 2 <consent> 3 {053 RATTLE failed to converge 
after 1000 iterations}
   while executing
"integrate 1000"
   ("while" body line 3)
   invoked from within
"while { $min < 1.50} {
   inter ljforcecap $cap
   integrate 1000
   set min [analyze mindist]
   incr cap 1000
puts -nonewline " $min "



#!/bin/sh
# tricking... the line after a these comments are interpreted as standard shell 
script \
    exec $ESPRESSO_SOURCE/Espresso $0 $*

puts "[code_info]"

#############################################################
#  System Parameter                                          #
#############################################################

set l_poly 44
set n_poly 240
set T 450
set lj1_epsilon [expr 0.09344/($T * 0.00198721)]    ; # Epsilon in KCal/mol 
divided by KbT in KCal/mol
set lj2_epsilon [expr 0.1455/($T * 0.00198721)]     ; # Epsilon in KCal/mol 
divided by KbT in KCal/mol
set lj3_epsilon [expr 0.2264/($T * 0.00198721)]     ; # Epsilon in KCal/mol 
divided by KbT in KCal/mol

set k1 [expr 120/($T * 0.00198721)]         ; # Stiffness/Bend in KCal/mol 
divided by KbT in KCal/mol
set k2 [expr 1.6/($T * 0.00198721)]         ; # Stiffness/Bend in KCal/mol 
divided by KbT in KCal/mol
set k3 [expr -0.867/($T * 0.00198721)]      ; # Stiffness/Bend in KCal/mol 
divided by KbT in KCal/mol
set k4 [expr 3.24/($T * 0.00198721)]        ; # Stiffness/Bend in KCal/mol 
divided by KbT in KCal/mol

#set sigma 4.5
set lj1_sig 4.5
set lj1_cut [expr 2*$lj1_sig]
set lj1_shift [expr 0.031]
set n_part  [expr $l_poly*$n_poly]
set box_length [expr 39.8886]
set bond_l [expr 1.53]
set shield 0.1

setmd box_l $box_length $box_length $box_length
setmd time_step 0.001
setmd skin 0.4

# Thermostat
integrate set nvt
set kT 1.0
thermostat langevin $kT 1

set vmd_script "vmd_start.script"

set vmd_online "no"
set vmd_offline "yes"
#############################################################
#  Interaction Setup                                          #
#############################################################

inter 1 rigid_bond $bond_l 1e-3 1e-3 
inter 2 angle $k1 [expr [PI]*0.61] ; # Cossquare Bond Angle Potential
inter 3 DIHEDRAL $k2 $k3 $k4           


inter 0 0 lj-gen $lj1_epsilon $lj1_sig $lj1_cut $lj1_shift 0 12 6 1 2 
inter 0 1 lj-gen $lj2_epsilon $lj1_sig $lj1_cut $lj1_shift 0 12 6 1 2 
inter 1 1 lj-gen $lj3_epsilon $lj1_sig $lj1_cut $lj1_shift 0 12 6 1 2 

puts "Interactions:\n   [inter]"

#############################################################
#  Polymer & Particle setup                                    #
#############################################################

polymer $n_poly $l_poly $bond_l mode SAW $shield bond 1

for { set i 0 } { $i < $n_part } { incr i } {
if {$i % 44 == 0} {part $i type 1} elseif {$i % 44 == 43} {part $i type 1} else 
{part $i type 0}
}


for {set i 0} { $i < $n_part -2} {incr i} {
if {$i % 44 != 42 && $i % 44 != 43} {part [expr $i+1] bond 2 [expr $i] [expr 
$i+2]}
}

for {set i 0} { $i < $n_part - 3} {incr i} {
if {$i % 44 != 41 && $i % 44 != 42 && $i % 44 != 43} {part [expr $i+1] bond 3 
[expr $i] [expr $i+2] [expr $i+3]}
}


set struct [open "|gzip -c - >phy.gz" w]
blockfile $struct write variable box_l
blockfile $struct write particles "id pos type " all
blockfile $struct write bonds all
close $struct


#############################################################
#  VMD Connection                                          #
#############################################################


writepsf "poll.psf" ; writepdb "poll.pdb"

if { $vmd_online == "yes" } then {
    for {set port 10000} { $port < 65000 } { incr port } {
        catch {imd connect $port} res
        if {$res == ""} break
    }
    if {$port == 65000} {
        error "could not open a port for IMD, please try offline"
    }
}

 if { $vmd_offline == "yes" || $vmd_online == "yes" } then {
    puts "Creating $vmd_script..."
    set vmdout_file [open $vmd_script "w"]
    if { $vmd_online == "yes" } then {
        puts $vmdout_file "mol load vtf poll_online.vtf"
        puts $vmdout_file "imd connect localhost $port"
    } elseif { $vmd_offline == "yes" } then {
        puts $vmdout_file "mol load vtf poll.vtf"
    }
      
    close $vmdout_file
    puts "$vmd_script finished."
}

# start vmd
if { $vmd_online == "yes" } then {
    set vtf_file [open "poll_online.vtf" w]
    writevsf $vtf_file
    puts $vtf_file 
    writevcf $vtf_file
    close $vtf_file

    # start VMD
    exec vmd -e $vmd_script &
    # wait for VMD to connect
    imd listen 100000
}


# create vtf-file
if { $vmd_offline == "yes" } then {
    set vtf_file [open "poll.vtf" w]
    writevsf $vtf_file
    # connect cup
    puts $vtf_file 
    writevcf $vtf_file
}



#############################################################
#  Warm up Integration                                     #
#############################################################

set min 0
set cap 1000
while { $min < 1.50} {
    inter ljforcecap $cap
    integrate 1000
    set min [analyze mindist]
    incr cap 1000
puts -nonewline " $min "

 if { $vmd_online=="yes" } { imd positions }
 if { $vmd_offline=="yes" } then { writevcf $vtf_file }
}
puts "Warmup finished. Minimal distance now $min"


#############################################################
#  Main Integration                                        #
#############################################################

set int_loop 20000000
set int_step 1


polyBlockWrite "poll.set"   {box_l time_step skin}
polyBlockWrite "poll.start" {time} {id pos type}

set obs [open "rg.dat" "w"]
setmd time 0; inter ljforcecap 0;set tmp_step 0
set i 0 
while { $i<$int_loop } {
set tmp_step [expr ($i+1)*$int_step]
set Temp [expr [analyze energy kinetic]/(($n_part*3 - 1935)/2.0)]
set p1 [analyze pressure total]
analyze set chains 0 $n_poly $l_poly

set rg [lindex [analyze rg] 0]
set re [lindex [analyze re] 2]
set Cn [expr $re /(43 * 1.53 * 1.53)] ;# Characteristic Ratio

if {$i % 100 == 0} {puts $obs "[setmd time] $rg  $re $Cn  $p1"}

 if { $i%1000==0 } {
      
        polyBlockWrite "poll[format %04d $i]" {time box_l time_step skin} {id 
pos type}
        
    }


set t [setmd time]
puts -nonewline "Step $tmp_step/[expr $int_step*$int_loop]  "
puts -nonewline  [format "    %.5f %.4f %.4f  %.4f  %.4f    %.4f...\r"  $t  $p1 
 $rg $re $Cn  $Temp]
integrate $int_step
if { $vmd == "yes" } { imd positions }
incr i

  if { $vmd_online=="yes" } { imd positions }
  if { $vmd_offline=="yes" } then { writevcf $vtf_file }
}
close $obs
if { $vmd_online=="yes" } { exec rm vmd_start.script }
if { $vmd_offline=="yes" } then { close $vtf_file }
puts "\nFinished"
exit