Enable work computation with thermostat heat
[namd.git] / lib / namdcph / namdcph / namdtcl.tcl
1 # namdtcl.tcl
2 #
3 # This file contains routines that mimic idiomatic NAMD behavior or otherwise
4 # re-implement utilities from the NAMD source code in pure Tcl.
5 #
6 # Check that this is being run through NAMD.
7 if {[catch numPes]} {
8     puts "Expected a NAMD-Tcl interpreter - exiting."
9     return 1
10 }
11
12 # =============================================================================
13 # Constants (These usually come from src/common.h)
14 # =============================================================================
15 variable BOLTZMANN 0.001987191 ;# in kcal/mol-K
16 variable PRESSUREFACTOR 69500.0 ;# convert kcal/mol to bar*A^3
17 variable PI 3.141592653589793
18 variable LN10 [expr {log(10)}]
19 variable LN2 [expr {log(2)}]
20
21 variable dummyXSC "# NAMD extended system configuration output file\n#\$LABELS step a_x a_y a_z b_x b_y b_z c_x c_y c_z o_x o_y o_z\n0 0 0 0 0 0 0 0 0 0 0 0 0\n"
22
23 # =============================================================================
24 # Standard energy callback
25 # =============================================================================
26 global EnergyLabels
27 global EnergyValues
28 global energyArray
29
30 # energyCallback
31 #
32 # Set the "EnergyLabels" and "EnergyValues" lists to contain the labels and 
33 # values for the most recent NAMD MD step, respectively. In practice this is 
34 # only meant to be used as an argument to the NAMD built-in "callback" 
35 # procedure.
36 #
37 # Arguments:
38 # ----------
39 # labels : list of strings
40 #   a list of energy term labels (e.g. {TS BOND ANGLE ... ELECT ...})
41 # values : list of floats
42 #   a list of energy term values corresponding to the above labels
43 #
44 # Returns:
45 # --------
46 # None
47 #
48 proc energyCallback {labels values} {
49     set ::EnergyLabels $labels
50     set ::EnergyValues $values
51     return
52 }
53
54 # storeEnergies
55 #
56 # Use the "EnergyLabels" and "EnergyValues" lists to populate the "energyArray"
57 # array. This is mostly a convenience function to cleanly and clearly retrieve
58 # energy information from NAMD. While NAMD will automatically call 
59 # energyCallback, this proc should be called each time "energyArray" is
60 # to be queried.
61 #
62 proc storeEnergies {} {
63     foreach label $::EnergyLabels value $::EnergyValues {
64         set ::energyArray($label) $value
65     }
66     return
67 }
68
69 # =============================================================================
70 # File handling
71 # =============================================================================
72 # Backup a file in the usual NAMD manner by appending ".BAK"
73 proc namdFileBackup {filename} {
74     if {[file exists $filename]} {
75         file copy -force $filename "$filename.BAK"
76     }
77     return
78 }
79
80 # =============================================================================
81 # Input checking (mimic functionality in src/SimParameters.C)
82 # =============================================================================
83 # Verify that the input value can be interpreted as a number. 
84 proc checkIsNumeric {name value} {
85     set IsNumeric 0
86     if {![catch {expr {abs($value)}}]} {
87         set IsNumeric 1
88     }
89     set value [string trimleft $value 0]
90     if {![catch {expr {abs($value)}}]} {
91         set IsNumeric 1
92     }
93     if {!$IsNumeric} {
94         abort "$name must be a numeric argument! (got $value)"
95     }
96     return 1
97 }
98
99 # Verify that the input value is a positive (nonzero) number.
100 proc checkIsPositive {name value} {
101     checkIsNumeric $name $value
102     if {$value <= 0} {
103         abort "$name must be greater than zero!"
104     }
105     return 1
106 }
107
108 # Verify that the input value is a non-negative (possibly zero) number.
109 proc checkIsNotNegative {name value} {
110     checkIsNumeric $name $value
111     if {$value < 0} {
112         abort "$name must be greater than or equal to zero!"
113     }
114     return 1
115 }
116
117 # Verify that a variable argument list is a multiple of the given number.
118 #
119 # Example: >> checkArglistIsMultiple argList 2
120 #          returns true if argList has an even number of elements.
121 #
122 proc checkArglistIsMultiple {argList multi} {
123     if {[expr {[llength $argList] % $multi}] != 0} {
124         abort "argument list must have multiple of $multi arguments!"
125     }
126     return 1
127 }
128
129 # =============================================================================
130 # Advanced parameter introspection
131 #
132 #   The following procs aim to solve the problem of keywords with redundant
133 # meaning. For example, multiple thermostats are available in NAMD and each has
134 # its own temperature variable. However, one generally does not care _how_ the
135 # temperature is regulated, but merely that it _is_ regulated and has a clear
136 # thermodynamic value. This is solved by having auxillary global tcl variables
137 # based on an exhaustive query of the keywords. For convenience, an extra 
138 # variable is also defined to give an explicit name to the thermostat so that
139 # consistent capitalization and formatting is achieved.
140 #
141 # Example: What is the temperature?
142 # >> # some MD operations
143 # >> getThermostat ;# populated the global variables
144 # >> set theCurrentTemperature [$::thermostatTempCmd] ;# query the temperature
145 # >> $::thermostatCmd off ;# turn off the termostat
146 #
147 # =============================================================================
148 global thermostatIsSet 0
149 global thermostatName ""
150 global thermostatCmd ""
151 global thermostatTempCmd ""
152
153 global barostatIsSet 0
154 global barostatName ""
155 global barostatCmd ""
156 global barostatPresCmd ""
157 global barostatTempCmd ""
158
159 # getThermostat
160 #
161 # This only requires that the thermostat itself has been invoked, not that a
162 # temperature has been set. By default, the much maligned Berendsen thermostat
163 # raises an error, since this is generally incompatible with canonical
164 # sampling.
165 #
166 # Return 1 if a thermostat is set, else return 0.
167 #
168 proc getThermostat {{forbidBerendsen true}} {
169     global thermostatIsSet
170     global thermostatName ""
171     global thermostatCmd ""
172     global thermostatTempCmd ""
173
174     # Only available in 2.12 development and 2.13+
175     if {![catch stochRescale]} {
176         if {[isset stochRescale] && [stochRescale]} {
177             set thermostatName "stochastic-rescaling"
178             set thermostatCmd stochRescale
179             set thermostatTempCmd stochRescaleTemp
180         }
181     }
182     if {[isset langevin] && [langevin]} {
183         set thermostatName "Langevin"
184         set thermostatCmd langevin
185         set thermostatTempCmd langevinTemp
186     } elseif {[isset loweAndersen] && [loweAndersen]} {
187         set thermostatName "Lowe-Andersen"
188         set thermostatCmd loweAndersen
189         set thermostatTempCmd loweAndersenTemp
190     } elseif {[expr {[reassignFreq] > 0}]} {
191         set thermostatName "Andersen (massive collisions)"
192         set thermostatCmd reassignFreq
193         set thermostatTempCmd reassignTemp
194     } elseif {[isset tCouple] && [tCouple]} {
195         if {$forbidBerendsen} {
196             abort "Berendsen thermostat does not work for canonical sampling."
197         } else {
198             set thermostatName "Berendsen"
199             set thermostatCmd tCouple
200             set thermostatTempCmd tCoupleTemp
201         }
202     } else {
203         set thermostatIsSet 0 
204         return 0
205     }
206     set thermostatIsSet 1
207     return 1 
208 }
209
210 # kBT
211 #
212 # Return kBT using the current thermostat temperature. If no thermostat is set
213 # then return 0.
214 #
215 proc kBT {} {
216     if {$::thermostatIsSet} {
217         return [expr {$::BOLTZMANN*[$::thermostatTempCmd]}]
218     }
219     return 0.0
220 }
221
222 # getBarostat
223 #
224 # This only requires that the barostat itself has been invoked, not that a
225 # pressure has been set. By default, the slightly less maligned Berendsen
226 # barostat raises an error, since this is generally incompatible with
227 # isobaric-isothermal sampling.
228 #
229 # Return 1 if a barostat is set, else return 0.
230 #
231 proc getBarostat {{forbidBerendsen true}} {
232     global barostatIsSet
233     global barostatName ""
234     global barostatCmd ""
235     global barostatPresCmd ""
236     global barostatTempCmd ""
237
238     if {[isset langevinPiston] && [langevinPiston]} {
239         set barostatName "Langevin Piston"
240         set barostatCmd langevinPiston
241         set barostatPresCmd langevinPistonTarget
242         set barostatTempCmd langevinPistonTemp
243     } elseif {[isset berendsenPressure] && [berendsenPressure]} {
244         if {$forbidBerendsen} {
245             abort "Berendsen barostat does not work for NpT sampling."
246         } else {
247             set barostatName "Berendsen"
248             set barostatCmd berendsenPressure
249             set barostatPresCmd berendsenPressureTarget
250         }
251     } else {
252         set barostatIsSet 0
253         return 0
254     }
255     set barostatIsSet 1
256     return 1 
257 }
258
259 # Alchemical Interactions
260
261 #   NAMD uses separate scaling schemes for bonded, electrostatic, and van der 
262 # Waals terms, each of which is determined by a separate parameter. This
263 # essentially permits staggered scaling so that some interactions can be turned
264 # off quickly (i.e. electrostatics) while others can be turned off more slowly
265 # (i.e. bonded and van der Waals). The following routines just re-implement 
266 # the NAMD internals so that only a single scaling parameter (alchLambda) needs
267 # to be tracked.
268 #
269 #   For convenience, all 6 scaling factors can also be obtained at the same 
270 # time - this is a bit opaque sometimes, but probably much cleaner than using
271 # separate variables for each energy component. 
272 #   There are 3 interaction types  (bonds, electrostatics, van der Waals - in 
273 # _that_ order) and two alchemical groups (1 and 2 - in _that_ order).
274 #
275 proc getLambdas {lambda1} {
276     set lambda2 [expr {1. - $lambda1}]
277     return [list [getBondLambda $lambda1] [getElecLambda $lambda1]\
278                  [getVdwLambda $lambda1] [getBondLambda $lambda2]\
279                  [getElecLambda $lambda2] [getVdwLambda $lambda2]\
280            ]
281 }
282
283 proc getElecLambda {lambda} {
284     set ls [expr [alchElecLambdaStart]]
285     return [expr {($lambda <= $ls) ? 0. : [expr {($lambda-$ls) / (1.-$ls)}]}]
286 }
287
288 proc getVdwLambda {lambda} {
289     set le [expr [alchVdwLambdaEnd]]
290     return [expr {($lambda >= $le) ? 1. : [expr {$lambda / $le}]}]
291 }
292
293 proc getBondLambda {lambda} {
294     set le [expr [alchBondLambdaEnd]]
295     return [expr {($lambda >= $le) ? 1. : [expr {$lambda / $le}]}]
296 }
297
298 # Return the alchemical energy components in the same list structure used by
299 # getLambdas. Note that energyArray should be passed _by name_, not by value.
300 #
301 # Example:
302 # callback energyCallback
303 # storeEnergies ;# this populates the global array energyArray
304 # set energyList [getAlchEnergies ::energyArray]
305
306 proc getAlchEnergies {energyArray} {
307     upvar 1 $energyArray energies
308     return [list $energies(BOND1) $energies(ELEC1) $energies(VDW1)\
309                  $energies(BOND2) $energies(ELEC2) $energies(VDW2)\
310            ]
311 }
312