/
el2d_nonneutral.tcl
140 lines (111 loc) · 3.95 KB
/
el2d_nonneutral.tcl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
# Copyright (C) 2010,2011,2012,2013,2014 The ESPResSo project
# Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
# Max-Planck-Institute for Polymer Research, Theory Group
#
# This file is part of ESPResSo.
#
# ESPResSo is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# ESPResSo is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
source "tests_common.tcl"
require_feature "ELECTROSTATICS"
require_feature "FFTW"
require_feature "PARTIAL_PERIODIC"
require_feature "CONSTRAINTS"
puts "-------------------------------------------"
puts "- Testcase el2d_nonneutral.tcl running on [format %02d [setmd n_nodes]] nodes: -"
puts "-------------------------------------------"
set epsilon 1e-3
thermostat off
setmd time_step 0.01
setmd skin 0.05
proc read_data {file} {
set f [open $file "r"]
while {![eof $f]} { blockfile $f read auto}
close $f
}
proc write_data {file} {
set f [open $file "w"]
blockfile $f write variable box_l
blockfile $f write tclvariable energy
blockfile $f write interactions
blockfile $f write particles {id pos q f}
close $f
}
if { [catch {
read_data "el2d_system.data"
for { set i 0 } { $i <= [setmd max_part] } { incr i } {
part $i q 1
}
set bl [setmd box_l]
constraint plate -height 0 -sigma [expr -[setmd max_part]/[lindex $bl 0]/[lindex $bl 1]]
# buffer the balanced node_grid for domain_decomposition/P3M
set balance_ng [setmd node_grid]
#########################################
#### check MMM2D first
#########################################
setmd periodic 1 1 0
cellsystem layered [expr 24/[setmd n_nodes]]
inter coulomb 2.0 mmm2d 1e-4
puts -nonewline "MMM2D node_grid=[setmd node_grid] "
flush stdout
integrate 0
for { set i 0 } { $i <= [setmd max_part] } { incr i } {
set F($i) [part $i pr f]
}
set energy [analyze energy total]
inter coulomb 0.0
puts "done"
if {[has_feature "FFTW"]} {
#########################################
#### check P3M + ELC
#########################################
cellsystem domain_decomposition
eval setmd node_grid $balance_ng
puts -nonewline "P3M+ELC node_grid=[setmd node_grid] "
flush stdout
setmd periodic 1 1 1
inter coulomb 2.0 p3m 16.3929265333 32 4 0.142069354322 9.78886014586e-05
inter coulomb epsilon metallic n_interpol 32768 mesh_off 0.5 0.5 0.5
inter coulomb elc 1e-4 [expr 0.1*[lindex [setmd box_l] 2]] -noneutralization
integrate 0
############## RMS force error for P3M + ELC
set rmsf 0
for { set i 0 } { $i <= [setmd max_part] } { incr i } {
set resF [part $i pr f]
set tgtF $F($i)
set dx [expr abs([lindex $resF 0] - [lindex $tgtF 0])]
set dy [expr abs([lindex $resF 1] - [lindex $tgtF 1])]
set dz [expr abs([lindex $resF 2] - [lindex $tgtF 2])]
set rmsf [expr $rmsf + $dx*$dx + $dy*$dy + $dz*$dz]
}
set rmsf [expr sqrt($rmsf/[setmd n_part])]
puts "rms force deviation $rmsf"
if { $rmsf > $epsilon } {
puts "P3M+ELC force error too large"
}
set cureng [lindex [analyze energy coulomb] 0]
set toteng [analyze energy total]
if { [expr abs($toteng - $cureng)] > $epsilon } {
error "system has unwanted energy contributions of [format %e [expr $toteng - $cureng]]"
}
set rel_eng_error [expr abs(($toteng - $energy)/$energy)]
puts "relative energy deviations: $rel_eng_error"
if { $rel_eng_error > $epsilon } {
error "relative energy error too large"
}
}
} res ] } {
error_exit $res
}
exit 0