Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Fetching contributors…

Cannot retrieve contributors at this time

file 122 lines (106 sloc) 3.88 kb
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
# Copyright (C) 2010,2011,2012 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/>.

#
#
# To generate the test system, use gen_harm.tcl.
#
source "tests_common.tcl"

require_feature "ADRESS" off

puts "------------------------------------------"
puts "- Testcase harm.tcl running on [format %02d [setmd n_nodes]] nodes: -"
puts "------------------------------------------"

setmd time_step 0.006
setmd skin 1.0
thermostat off

set epsilon 1e-4

proc read_data {file} {
    if { [string compare [lindex [split $file "."] end] "gz"]==0 } { set f [open "|gzip -cd $file" r]
    } else { set f [open "$file" "r"] }
    while {![eof $f]} { blockfile $f read auto}
   if { [catch { close $f } fid] } { puts "Error while closing $file caught: $fid." }
}

if { [catch {
    read_data "harm_system.data.gz"

    for { set i 0 } { $i <= [setmd max_part] } { incr i } {
set F($i) [part $i pr f]
    }
    # to ensure force recalculation
    invalidate_system
    integrate 0
    # here you can create the necessary snapshot
    # write_data "harm_system.data2"

    # ensures that no other forces are on
    set cureng [analyze energy bonded 0 0]
    # tbrs
    set curprs [lindex [analyze pressure bonded 0 0] 0]
    ############## end

    set toteng [analyze energy total]
    set totprs [analyze pressure total]

    if { [expr abs($toteng - $cureng)] > $epsilon } {
error "system has unwanted energy contributions"
    }
    if { [expr abs($totprs - $curprs)] > $epsilon } {
error "system has unwanted pressure contributions"
    }

    set rel_eng_error [expr abs(($toteng - $energy)/$energy)]
    puts "relative energy deviations: $rel_eng_error ($toteng / $energy)"
    if { $rel_eng_error > $epsilon } {
error "relative energy error too large"
    }

    set rel_prs_error [expr abs(($totprs - $pressure)/$pressure)]
    puts "relative pressure deviations: $rel_prs_error ($totprs / $pressure)"
    if { $rel_prs_error > $epsilon } {
error "relative pressure error too large"
    }

    set maxdx 0
    set maxpx 0
    set maxdy 0
    set maxpy 0
    set maxdz 0
    set maxpz 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])]

if { $dx > $maxdx} {
set maxdx $dx
set maxpx $i
}
if { $dy > $maxdy} {
set maxdy $dy
set maxpy $i
}
if { $dz > $maxdz} {
set maxdz $dz
set maxpz $i
}
    }
    puts "maximal force deviation in x $maxdx for particle $maxpx, in y $maxdy for particle $maxpy, in z $maxdz for particle $maxpz"
    if { $maxdx > $epsilon || $maxdy > $epsilon || $maxdz > $epsilon } {
if { $maxdx > $epsilon} {puts "force of particle $maxpx: [part $maxpx pr f] != $F($maxpx)"}
if { $maxdy > $epsilon} {puts "force of particle $maxpy: [part $maxpy pr f] != $F($maxpy)"}
if { $maxdz > $epsilon} {puts "force of particle $maxpz: [part $maxpz pr f] != $F($maxpz)"}
error "force error too large"
    }
} res ] } {
    error_exit $res
}

exit 0
Something went wrong with that request. Please try again.