|
(You
may want to widen your browser's window until the
lines of text in the boxes are no longer
wrapped.)
Below
is a simple example. The script starts with the
required importation of the
yammp
module and call to the
setup()
function. Following that, the nonbond functions,
the Monte Carlo and the Molecular Dynamics methods
are modified. Next, the Monte Carlo method is
called and the return value is captured in a
variable named discard which is then ignored.
Finally, a complex set of operations are carried
out two times. The set starts by setting the
starting temperature to 1 K, the final temperature
to 250 K, generating new velocities at the starting
temperature and running 10000 steps of Molecular
Dynamics. This is followed by setting the starting
and final temperatures to 250 K, scaling the
prevailing velocities to the starting temperature
and running 25000 steps of Molecular Dynamics. The
set continues by setting the final temperature to 1
K, scaling the velocities to the starting
temperature of 250 K and running 500000 steps of
Molecular Dynamics. At the end of the set, the
coordinates are written to the output archive
file.
from Taro.Yammp import yammp
discard = yammp.setup( descriptor = "tRNA.ds", archive = "tRNA0.ar", outarchive = "tRNA.ar" )
yammp.editnonbonds( vlatscale = 5 )
yammp.editmontecarlo( report = 100 )
yammp.editdynamics(
step_size = 0.001,
report = 25000,
scale_velocities = 25,
scale_type = "berendsen",
thermal = 0.01 )
discard = yammp.montecarlo( 1000 )
for i in range( 2 ):
yammp.setenviron( start_temperature = 1, final_temperature = 250, set_velocities = "new" )
discard = yammp.dynamics( 10000 )
yammp.setenviron( temperature = 250, set_velocities = "scale" )
discard = yammp.dynamics( 25000 )
yammp.setenviron( start_temperature = 250, final_temperature = 1, set_velocities = "scale" )
discard = yammp.dynamics( 500000 )
discard = yammp.archive()
|
In
the following example, we write a module named
anneal
and the module contains a function named
withall().
The file is named anneal.py so that Python can
recognize it as a module file.
The
top of the file contains a string of three lines
enclosed by a pair of triple quotation marks. This
is called a module documentation string and it
serves to document the code and it is accessible to
the user of the module. We could put the note in a
comment (starting with #) but comments are not
accessible except by reading the source code. The
documentation string refer to a second function
withmc()
which we have removed to reduce clutter.
Following
the module documentation string we import the
yammp
module from the Taro package.
Next,
we define the function
withall().
Note the triple quoted string. This is the function
documentation string and it serves a similar
purpose to the the module documentation
string.
We
name the formal arguments in the function header.
On the last three arguments we set default values.
This means that this function can be called with as
few as one argument (provided the default values
for the other arguments are acceptable).
After
setting some internal variables (we name them with
a leading underscore), we check the values of the
argument. For example, we make sure that the final
temperature is lower than the initial, swapping
them if not. This is probably not a good idea but
this is just an example. We then take the number of
steps that the user has supplied and portion it up
unequally among three molecular mechanics
methods.
In
the first portion of steps, we actually use half
our budget in Monte Carlo, heating the molecule at
the starting temperature. Monte Carlo is a good
method for this purpose because with Molecular
Dynamics we would have to worry about stability and
would have to control the heating
process.
We
than use most of our budget of steps in Molecular
Dynamics in several portions. In each portion we
bring the temperature down midway between the
initial and the current target final temperature,
i.e, we are halving the temperature difference.
This gives is a rapid drop in temperature at the
start and the temperature tails off slowly. Again,
this is contrived to make the example interesting.
Note that we cannot reach the target final
temperature specified by the user.
In
the final portion we again use only half our budget
of steps in Energy Minimization. At this stage, the
molecule should be at a very low energy and we
really need only a few Minimization steps to polish
the structure.
anneal.py
|
"""anneal: this module contains a collection of simulated annealing functions
withmc() -- simple procedure using only Monte Carlo
withall() -- a more complex procedure using MC, MD and EM"""
from Taro.Yammp import yammp
def withall( steps, start=500.0, final=0.1, stepsize=0.001 ):
"""withall( steps, start=500, final=0.1, stepsize=0.001 ) -- annealing in
several parts: first hold at starting temperature using Monte Carlo
then switching to Molecular Dynamics to follow a decay temperature curve
from 'start' to 'finish', finished with Energy Minimization"""
# ----- internal settings
_diminish = 0.5
_minsteps = 10000
_numparts = 10
# ----- check supplied values
if start < final:
print "## start temperature (%.8g K) lower than final (%.8g K), will swap\n"\
% ( start, final )
start, final = final, start
if steps < _minsteps:
print "## number of steps (%d) is too few, will use %d\n"\
% ( steps, _minsteps )
steps = _minsteps
_partsteps = steps / _numparts
_partprint = _partsteps / 2
# ----- print the energy and rms gradient at the start
E, RMS_grad = yammp.energy_rmsgrad()
print "\n>>> Simulated Annealing using MC, MD & EM, Initial Energy",
print "%.6g kcal/mol, RMS Gradient %.6g kcal/mol/A\n\n" % ( E, RMS_grad )
# ----- in the first part, heat to the starting temperature
print "<< First Part: Holding at %.8g K over %d cycles\n" % ( start, _partsteps / 2 )
yammp.setenviron( temperature = start )
yammp.editmontecarlo(
report = _partprint / 2,
update_list = 100 )
yammp.montecarlo( _partsteps / 2 )
# ----- diminish the temperature drop in each of the remaining parts
yammp.editdynamics(
step_size = stepsize,
report = _partprint,
scale_velocities = 25,
scale_type = "berendsen",
thermal = 10 * stepsize )
yammp.setenviron(
temperature = start,
set_velocities = "new" )
delta = ( start - final ) * _diminish
thisstart = start
for part in range( _numparts - 2 ):
thisfinal = thisstart - delta
print "<< Part %d: Cooling from %.8g K to %.8g K over %d steps\n"\
% ( part + 2, thisstart, thisfinal, _partsteps )
yammp.setenviron(
start_temperature = thisstart,
final_temperature = thisfinal,
set_velocities = "scale" )
yammp.dynamics( _partsteps )
thisstart = thisfinal
delta = ( thisstart - final ) * _diminish
# ----- polish off with energy minimization
print "<< Final Part: Minimize to a maximum of %d steps\n" % ( _partsteps / 2 )
yammp.editlinesearch( iterations = 3 )
yammp.editminimizers(
report = _partprint,
update_list = 1000 )
yammp.minimize( _partsteps / 2 )
E, RMS_grad = yammp.energy_rmsgrad()
print "\n>>> Simulated Annealing completed, Final Energy",
print "%.6g kcal/mol, RMS Gradient %.6g kcal/mol/A\n" % ( E, RMS_grad )
|
If
we have a complex sequence of calculations that we
use repeatedly we can package the calculations into
a function. If the sequence can be used with
different settings, we have an even more useful
function. As we develop more functions we can
package related functions into modules.
Finally,
this is one way of using the module we have just
developed. We first import the
yammp
module since we need to use the
setup()
and archive()
functions. Next, we import the
anneal
module.
We
have to call setup()
before anything else. Then we change a setting in
the nonbond term of the potential energy function.
Finally, a for loop to carry out a sequence of
actions twice. In each sequence we call the
withall()
function and on return from the function we save
the coordinates to an output archive file. Note
that we accept the default step size in the
withall()
function.
from Taro.Yammp import yammp
import anneal
yammp.setup(
descriptor = "dna.ds",
archive = "dna0.ar",
outarchive = "dna1.ar" )
yammp.editnonbonds( vlatscale = 5 )
for i in range( 2 ):
anneal.withall(
steps = 500000,
start = 400,
final = 0 )
yammp.archive()
|
More
examples?
|