Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78663362
temper.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Aug 22, 06:42
Size
11 KB
Mime Type
text/x-c
Expires
Sat, Aug 24, 06:42 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20082838
Attached To
rLAMMPS lammps
temper.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mark Sears (SNL)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "temper.h"
#include "universe.h"
#include "domain.h"
#include "atom.h"
#include "update.h"
#include "integrate.h"
#include "modify.h"
#include "compute.h"
#include "force.h"
#include "output.h"
#include "thermo.h"
#include "fix.h"
#include "random_park.h"
#include "finish.h"
#include "timer.h"
#include "memory.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
// #define TEMPER_DEBUG 1
/* ---------------------------------------------------------------------- */
Temper
::
Temper
(
LAMMPS
*
lmp
)
:
Pointers
(
lmp
)
{}
/* ---------------------------------------------------------------------- */
Temper
::~
Temper
()
{
MPI_Comm_free
(
&
roots
);
if
(
ranswap
)
delete
ranswap
;
delete
ranboltz
;
delete
[]
set_temp
;
delete
[]
temp2world
;
delete
[]
world2temp
;
delete
[]
world2root
;
}
/* ----------------------------------------------------------------------
perform tempering with inter-world swaps
------------------------------------------------------------------------- */
void
Temper
::
command
(
int
narg
,
char
**
arg
)
{
if
(
universe
->
nworlds
==
1
)
error
->
all
(
FLERR
,
"Must have more than one processor partition to temper"
);
if
(
domain
->
box_exist
==
0
)
error
->
all
(
FLERR
,
"Temper command before simulation box is defined"
);
if
(
narg
!=
6
&&
narg
!=
7
)
error
->
universe_all
(
FLERR
,
"Illegal temper command"
);
int
nsteps
=
force
->
inumeric
(
FLERR
,
arg
[
0
]);
nevery
=
force
->
inumeric
(
FLERR
,
arg
[
1
]);
double
temp
=
force
->
numeric
(
FLERR
,
arg
[
2
]);
for
(
whichfix
=
0
;
whichfix
<
modify
->
nfix
;
whichfix
++
)
if
(
strcmp
(
arg
[
3
],
modify
->
fix
[
whichfix
]
->
id
)
==
0
)
break
;
if
(
whichfix
==
modify
->
nfix
)
error
->
universe_all
(
FLERR
,
"Tempering fix ID is not defined"
);
seed_swap
=
force
->
inumeric
(
FLERR
,
arg
[
4
]);
seed_boltz
=
force
->
inumeric
(
FLERR
,
arg
[
5
]);
my_set_temp
=
universe
->
iworld
;
if
(
narg
==
7
)
my_set_temp
=
force
->
inumeric
(
FLERR
,
arg
[
6
]);
// swap frequency must evenly divide total # of timesteps
if
(
nevery
==
0
)
error
->
universe_all
(
FLERR
,
"Invalid frequency in temper command"
);
nswaps
=
nsteps
/
nevery
;
if
(
nswaps
*
nevery
!=
nsteps
)
error
->
universe_all
(
FLERR
,
"Non integer # of swaps in temper command"
);
// fix style must be appropriate for temperature control
if
((
strcmp
(
modify
->
fix
[
whichfix
]
->
style
,
"nvt"
)
!=
0
)
&&
(
strcmp
(
modify
->
fix
[
whichfix
]
->
style
,
"langevin"
)
!=
0
)
&&
(
strcmp
(
modify
->
fix
[
whichfix
]
->
style
,
"temp/berendsen"
)
!=
0
)
&&
(
strcmp
(
modify
->
fix
[
whichfix
]
->
style
,
"temp/rescale"
)
!=
0
))
error
->
universe_all
(
FLERR
,
"Tempering temperature fix is not valid"
);
// setup for long tempering run
update
->
whichflag
=
1
;
update
->
nsteps
=
nsteps
;
update
->
beginstep
=
update
->
firststep
=
update
->
ntimestep
;
update
->
endstep
=
update
->
laststep
=
update
->
firststep
+
nsteps
;
if
(
update
->
laststep
<
0
)
error
->
all
(
FLERR
,
"Too many timesteps"
);
lmp
->
init
();
// local storage
me_universe
=
universe
->
me
;
MPI_Comm_rank
(
world
,
&
me
);
nworlds
=
universe
->
nworlds
;
iworld
=
universe
->
iworld
;
boltz
=
force
->
boltz
;
// pe_compute = ptr to thermo_pe compute
// notify compute it will be called at first swap
int
id
=
modify
->
find_compute
(
"thermo_pe"
);
if
(
id
<
0
)
error
->
all
(
FLERR
,
"Tempering could not find thermo_pe compute"
);
Compute
*
pe_compute
=
modify
->
compute
[
id
];
pe_compute
->
addstep
(
update
->
ntimestep
+
nevery
);
// create MPI communicator for root proc from each world
int
color
;
if
(
me
==
0
)
color
=
0
;
else
color
=
1
;
MPI_Comm_split
(
universe
->
uworld
,
color
,
0
,
&
roots
);
// RNGs for swaps and Boltzmann test
// warm up Boltzmann RNG
if
(
seed_swap
)
ranswap
=
new
RanPark
(
lmp
,
seed_swap
);
else
ranswap
=
NULL
;
ranboltz
=
new
RanPark
(
lmp
,
seed_boltz
+
me_universe
);
for
(
int
i
=
0
;
i
<
100
;
i
++
)
ranboltz
->
uniform
();
// world2root[i] = global proc that is root proc of world i
world2root
=
new
int
[
nworlds
];
if
(
me
==
0
)
MPI_Allgather
(
&
me_universe
,
1
,
MPI_INT
,
world2root
,
1
,
MPI_INT
,
roots
);
MPI_Bcast
(
world2root
,
nworlds
,
MPI_INT
,
0
,
world
);
// create static list of set temperatures
// allgather tempering arg "temp" across root procs
// bcast from each root to other procs in world
set_temp
=
new
double
[
nworlds
];
if
(
me
==
0
)
MPI_Allgather
(
&
temp
,
1
,
MPI_DOUBLE
,
set_temp
,
1
,
MPI_DOUBLE
,
roots
);
MPI_Bcast
(
set_temp
,
nworlds
,
MPI_DOUBLE
,
0
,
world
);
// create world2temp only on root procs from my_set_temp
// create temp2world on root procs from world2temp,
// then bcast to all procs within world
world2temp
=
new
int
[
nworlds
];
temp2world
=
new
int
[
nworlds
];
if
(
me
==
0
)
{
MPI_Allgather
(
&
my_set_temp
,
1
,
MPI_INT
,
world2temp
,
1
,
MPI_INT
,
roots
);
for
(
int
i
=
0
;
i
<
nworlds
;
i
++
)
temp2world
[
world2temp
[
i
]]
=
i
;
}
MPI_Bcast
(
temp2world
,
nworlds
,
MPI_INT
,
0
,
world
);
// if restarting tempering, reset temp target of Fix to current my_set_temp
if
(
narg
==
7
)
{
double
new_temp
=
set_temp
[
my_set_temp
];
modify
->
fix
[
whichfix
]
->
reset_target
(
new_temp
);
}
// setup tempering runs
int
i
,
which
,
partner
,
swap
,
partner_set_temp
,
partner_world
;
double
pe
,
pe_partner
,
boltz_factor
,
new_temp
;
if
(
me_universe
==
0
&&
universe
->
uscreen
)
fprintf
(
universe
->
uscreen
,
"Setting up tempering ...
\n
"
);
update
->
integrate
->
setup
();
if
(
me_universe
==
0
)
{
if
(
universe
->
uscreen
)
{
fprintf
(
universe
->
uscreen
,
"Step"
);
for
(
int
i
=
0
;
i
<
nworlds
;
i
++
)
fprintf
(
universe
->
uscreen
,
" T%d"
,
i
);
fprintf
(
universe
->
uscreen
,
"
\n
"
);
}
if
(
universe
->
ulogfile
)
{
fprintf
(
universe
->
ulogfile
,
"Step"
);
for
(
int
i
=
0
;
i
<
nworlds
;
i
++
)
fprintf
(
universe
->
ulogfile
,
" T%d"
,
i
);
fprintf
(
universe
->
ulogfile
,
"
\n
"
);
}
print_status
();
}
timer
->
init
();
timer
->
barrier_start
();
for
(
int
iswap
=
0
;
iswap
<
nswaps
;
iswap
++
)
{
// run for nevery timesteps
update
->
integrate
->
run
(
nevery
);
// compute PE
// notify compute it will be called at next swap
pe
=
pe_compute
->
compute_scalar
();
pe_compute
->
addstep
(
update
->
ntimestep
+
nevery
);
// which = which of 2 kinds of swaps to do (0,1)
if
(
!
ranswap
)
which
=
iswap
%
2
;
else
if
(
ranswap
->
uniform
()
<
0.5
)
which
=
0
;
else
which
=
1
;
// partner_set_temp = which set temp I am partnering with for this swap
if
(
which
==
0
)
{
if
(
my_set_temp
%
2
==
0
)
partner_set_temp
=
my_set_temp
+
1
;
else
partner_set_temp
=
my_set_temp
-
1
;
}
else
{
if
(
my_set_temp
%
2
==
1
)
partner_set_temp
=
my_set_temp
+
1
;
else
partner_set_temp
=
my_set_temp
-
1
;
}
// partner = proc ID to swap with
// if partner = -1, then I am not a proc that swaps
partner
=
-
1
;
if
(
me
==
0
&&
partner_set_temp
>=
0
&&
partner_set_temp
<
nworlds
)
{
partner_world
=
temp2world
[
partner_set_temp
];
partner
=
world2root
[
partner_world
];
}
// swap with a partner, only root procs in each world participate
// hi proc sends PE to low proc
// lo proc make Boltzmann decision on whether to swap
// lo proc communicates decision back to hi proc
swap
=
0
;
if
(
partner
!=
-
1
)
{
if
(
me_universe
>
partner
)
MPI_Send
(
&
pe
,
1
,
MPI_DOUBLE
,
partner
,
0
,
universe
->
uworld
);
else
MPI_Recv
(
&
pe_partner
,
1
,
MPI_DOUBLE
,
partner
,
0
,
universe
->
uworld
,
MPI_STATUS_IGNORE
);
if
(
me_universe
<
partner
)
{
boltz_factor
=
(
pe
-
pe_partner
)
*
(
1.0
/
(
boltz
*
set_temp
[
my_set_temp
])
-
1.0
/
(
boltz
*
set_temp
[
partner_set_temp
]));
if
(
boltz_factor
>=
0.0
)
swap
=
1
;
else
if
(
ranboltz
->
uniform
()
<
exp
(
boltz_factor
))
swap
=
1
;
}
if
(
me_universe
<
partner
)
MPI_Send
(
&
swap
,
1
,
MPI_INT
,
partner
,
0
,
universe
->
uworld
);
else
MPI_Recv
(
&
swap
,
1
,
MPI_INT
,
partner
,
0
,
universe
->
uworld
,
MPI_STATUS_IGNORE
);
#ifdef TEMPER_DEBUG
if
(
me_universe
<
partner
)
printf
(
"SWAP %d & %d: yes = %d,Ts = %d %d, PEs = %g %g, Bz = %g %g
\n
"
,
me_universe
,
partner
,
swap
,
my_set_temp
,
partner_set_temp
,
pe
,
pe_partner
,
boltz_factor
,
exp
(
boltz_factor
));
#endif
}
// bcast swap result to other procs in my world
MPI_Bcast
(
&
swap
,
1
,
MPI_INT
,
0
,
world
);
// rescale kinetic energy via velocities if move is accepted
if
(
swap
)
scale_velocities
(
partner_set_temp
,
my_set_temp
);
// if my world swapped, all procs in world reset temp target of Fix
if
(
swap
)
{
new_temp
=
set_temp
[
partner_set_temp
];
modify
->
fix
[
whichfix
]
->
reset_target
(
new_temp
);
}
// update my_set_temp and temp2world on every proc
// root procs update their value if swap took place
// allgather across root procs
// bcast within my world
if
(
swap
)
my_set_temp
=
partner_set_temp
;
if
(
me
==
0
)
{
MPI_Allgather
(
&
my_set_temp
,
1
,
MPI_INT
,
world2temp
,
1
,
MPI_INT
,
roots
);
for
(
i
=
0
;
i
<
nworlds
;
i
++
)
temp2world
[
world2temp
[
i
]]
=
i
;
}
MPI_Bcast
(
temp2world
,
nworlds
,
MPI_INT
,
0
,
world
);
// print out current swap status
if
(
me_universe
==
0
)
print_status
();
}
timer
->
barrier_stop
();
update
->
integrate
->
cleanup
();
Finish
finish
(
lmp
);
finish
.
end
(
1
);
update
->
whichflag
=
0
;
update
->
firststep
=
update
->
laststep
=
0
;
update
->
beginstep
=
update
->
endstep
=
0
;
}
/* ----------------------------------------------------------------------
scale kinetic energy via velocities a la Sugita
------------------------------------------------------------------------- */
void
Temper
::
scale_velocities
(
int
t_partner
,
int
t_me
)
{
double
sfactor
=
sqrt
(
set_temp
[
t_partner
]
/
set_temp
[
t_me
]);
double
**
v
=
atom
->
v
;
int
nlocal
=
atom
->
nlocal
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
v
[
i
][
0
]
=
v
[
i
][
0
]
*
sfactor
;
v
[
i
][
1
]
=
v
[
i
][
1
]
*
sfactor
;
v
[
i
][
2
]
=
v
[
i
][
2
]
*
sfactor
;
}
}
/* ----------------------------------------------------------------------
proc 0 prints current tempering status
------------------------------------------------------------------------- */
void
Temper
::
print_status
()
{
if
(
universe
->
uscreen
)
{
fprintf
(
universe
->
uscreen
,
BIGINT_FORMAT
,
update
->
ntimestep
);
for
(
int
i
=
0
;
i
<
nworlds
;
i
++
)
fprintf
(
universe
->
uscreen
,
" %d"
,
world2temp
[
i
]);
fprintf
(
universe
->
uscreen
,
"
\n
"
);
}
if
(
universe
->
ulogfile
)
{
fprintf
(
universe
->
ulogfile
,
BIGINT_FORMAT
,
update
->
ntimestep
);
for
(
int
i
=
0
;
i
<
nworlds
;
i
++
)
fprintf
(
universe
->
ulogfile
,
" %d"
,
world2temp
[
i
]);
fprintf
(
universe
->
ulogfile
,
"
\n
"
);
fflush
(
universe
->
ulogfile
);
}
}
Event Timeline
Log In to Comment