Re: [ESPResSo-users] WARNING: pressure calculated, but MMM2D pressure not implemented

[Axel Arnold]

Hi Uday and Georg,

long ago we were discussing about implementing pressure for MMM2D and/or ELC. First of all, the simple dU/dV approach only works for 3d systems and only for the isotropic pressure, so one actually needs to compute the full pressure tensor.

In principle one can derive that from the MMM2D formulas, but practically that is quite tedious: at the heart of MMM2D/ELC there are trigonometric addition theorems. When deriving the cross terms of the pressure tensor, you get so many contributions that we found it impossible to implement the algorithm by hand, and writing a program to generate the code seemed somewhat overdoing it. Things got even more difficult because the current MMM2D/ELC also handles dielectric contrasts, which made the formulas even more awkward...

Best,

Axel

Am 26.09.2016 um 19:10 schrieb Georg Rempfer <address@hidden>:

Hello Padidela,

this just means that when you run 'analyze pressure total', the resulting value does not contain the pressure contribution due to the electrostatic interactions. There is nothing wrong with the dynamics of the system, this only influences this calculated observable.

Do you need to know the pressure? I assume it's not a problem in principle to include this contribution in the calculation, but so far no one has implemented the necessary code.

Greetings,

Georg

On Fri, Sep 16, 2016 at 3:39 PM, padidela uday <address@hidden> wrote:

Dear all,

I am trying to implement charged wall constraints in one of my scripts along with the particles between the constraints. But however I land up in few troubles with it.

After equilibrating the system a warning is displayed on screen during sampling

calculating pressure for electrostatics method that does have it implemented

WARNING: pressure calculated, but MMM2D pressure not implemented

Is there something wrong with the code could you help me to handle this?

Thank you for your time

set n_solute 200
set n_solvent 200
set top_wall 100
set bottom_wall 100
set n_wall 200
set n_part [expr ($n_solute + $n_solvent + $top_wall + $bottom_wall)]
set box_l 200
set sig1 2.0
set box_lxy [expr sqrt(2) * $box_l]
set box_lz [expr 0.5*$box_l]

setmd box_l $box_lxy $box_lxy [expr $box_lz +  $sig1]
setmd periodic 1 1 0 

for {set i 0} { $i < 200 } {incr i} {
    set posx [expr $box_lxy*[t_random]]
    set posy [expr $box_lxy*[t_random]]
    set posz [expr $box_lz*[t_random] + $sig1]
    set vx [gauss_random]
    set vy [gauss_random]
    set vz [gauss_random]
    part $i pos $posx $posy $posz type 0 v $vx $vy $vz q -2 mass 1.0
}

for {set i 200} { $i < 400 } {incr i} {
    set posx [expr $box_lxy*[t_random]]
    set posy [expr $box_lxy*[t_random]]
    set posz [expr $box_lz*[t_random] + $sig1]
    set vx [gauss_random]
    set vy [gauss_random]
    set vz [gauss_random]
    part $i pos $posx $posy $posz type 1 v $vx $vy $vz q 1 mass 1.0
}

#wall at top
constraint wall normal 0 0  1 dist $sig1 type 2

for {set i 400} { $i < 500 } {incr i} {
    set posx [expr $box_lxy*[t_random]]
    set posy [expr $box_lxy*[t_random]]
    set posz [expr ($box_lz + $sig1)]

    part  [expr $i + $one_wall] pos $posx $posy $posz type 2  q 1
}

constraint wall normal 0 0 -1 dist [expr -$box_lz - $box_lz] type 2

#wall at bottom

for {set i 500} { $i < 600 } {incr i} {
    set posx [expr $box_lxy*[t_random]]
    set posy [expr $box_lxy*[t_random]]
    set posz [expr ($box_lz - $box_lz)*[t_random]]
    part  [expr $i + $two_wall] pos $posx $posy $posz type 2  q 1
}

set sigma [expr -0.25*$n_wall/($box_lxy*$box_lxy)]
constraint plate height 0 sigma $sigma
constraint plate height [expr $box_lz + $sig1] sigma $sigma

setmd time_step 0.00002624
setmd skin 0.3
set temp 4.9576
set gamma 1
thermostat langevin $temp $gamma

#  POTENTIAL
set sig1 2.0
set sig2 1.15
set sig3 0.3
set eps1 1.0
set eps2 0.7
set eps3 0.5
set cut1 [expr 1.12246*$sig1]
set cut2 [expr 1.12246*$sig2]
set cut3 [expr 1.12246*$sig3]
set shift1 [expr 0.25*$eps3]
set shift2 [expr 0.25*$eps3]
set shift3 [expr 0.25*$eps3]
inter 0 0 lennard-jones $eps3 $sig1 $cut1 $shift1  0 
inter 1 0 lennard-jones $eps3 $sig2 $cut2 $shift2  0 
inter 1 1 lennard-jones $eps3 $sig3 $cut3 $shift3  0 
inter 2 0 lennard-jones $eps3 $sig1 $cut1 $shift1  0
inter 2 1 lennard-jones $eps3 $sig2 $cut2 $shift2  0
#inter 2 2 lennard-jones $eps3 $sig3 $cut3 $shift3  0 

cellsystem layered 3
inter coulomb 1.0 mmm2d 1e-4

if { [regexp "ROTATION" [code_info]] } {
      set deg_free 6
} {
      set deg_free 3
}

# prepare vmd connection
if { $vmd_output=="yes" } {
   prepare_vmd_connection "vmd"
    imd listen 10000
}

inter ljforcecap individual
set inter_steps 10000
for {set cap 10} {$cap < 100} {incr  10} {
       set rad [expr 1.0 - 0.5*$cap/10.0]
       set lb [expr 1.0 * $cap / 10.0]
       inter 0 0 lennard-jones $eps3 $sig1 $cut1 $shift1  0  $rad
       inter 1 0 lennard-jones $eps3 $sig2 $cut2 $shift2  0  $rad
       inter 1 1 lennard-jones $eps3 $sig3 $cut3 $shift3  0  $rad
       inter coloumb $lb mmm2d 1e-4
       set temp [expr [analyze energy kinetic]/(($deg_free/2.0)*$n_part)]
    puts "t=[setmd time] E=[analyze energy total] T=$temp"
         if {$vmd_output=="yes"} {imd positions}
    inter ljforcecap $cap;
    integrate $inter_steps
}
inter ljforcecap 0
inter coulomb 1.0 mmm2d 1e-4

set h [open "production.data" "w"]
set g [open "analysis.data" "w"]
set n_part [expr ($n_solute + $n_solvent + $top_wall + $bottom_wall)]
for {set i 0} { $i <4000 } {incr i} {
    set temp [expr [analyze energy kinetic]/(($deg_free/2.0)*$n_part)]
        puts $h " $i [setmd time] [analyze energy kinetic] $temp"
        integrate 200

    set f [open "config_$i" "w"]
    blockfile $f write tclvariable {box_lxy box_lz}
    blockfile $f write variable box_l
    blockfile $f write particles {id pos type v f}
    set temp [expr [analyze energy kinetic]/(($deg_free/2.0)*$n_part)]
    puts $f " \ { energy [analyze energy total] $temp}"
    puts $g " [analyze energy kinetic] [analyze pressure total] [analyze energy total] $temp"

    close $f
}
close $g

With Kind Regards

Uday