Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80986726
colvarproxy_lammps.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
Tue, Sep 3, 15:20
Size
14 KB
Mime Type
text/x-c
Expires
Thu, Sep 5, 15:20 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
20493673
Attached To
rLAMMPS lammps
colvarproxy_lammps.cpp
View Options
#include <mpi.h>
#include "lammps.h"
#include "atom.h"
#include "error.h"
#include "output.h"
#include "random_park.h"
#include "fix_colvars.h"
#include "colvarmodule.h"
#include "colvaratoms.h"
#include "colvarproxy.h"
#include "colvarproxy_lammps.h"
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <cerrno>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <sstream>
#include <string>
#define HASH_FAIL -1
////////////////////////////////////////////////////////////////////////
// local helper functions
// safely move filename to filename.extension
static
int
my_backup_file
(
const
char
*
filename
,
const
char
*
extension
)
{
struct
stat
sbuf
;
if
(
stat
(
filename
,
&
sbuf
)
==
0
)
{
if
(
!
extension
)
extension
=
".BAK"
;
char
*
backup
=
new
char
[
strlen
(
filename
)
+
strlen
(
extension
)
+
1
];
strcpy
(
backup
,
filename
);
strcat
(
backup
,
extension
);
#if defined(_WIN32) && !defined(__CYGWIN__)
remove
(
backup
);
#endif
if
(
rename
(
filename
,
backup
))
{
char
*
sys_err_msg
=
strerror
(
errno
);
if
(
!
sys_err_msg
)
sys_err_msg
=
(
char
*
)
"(unknown error)"
;
fprintf
(
stderr
,
"Error renaming file %s to %s: %s
\n
"
,
filename
,
backup
,
sys_err_msg
);
delete
[]
backup
;
return
COLVARS_ERROR
;
}
delete
[]
backup
;
}
return
COLVARS_OK
;
}
////////////////////////////////////////////////////////////////////////
colvarproxy_lammps
::
colvarproxy_lammps
(
LAMMPS_NS
::
LAMMPS
*
lmp
,
const
char
*
inp_name
,
const
char
*
out_name
,
const
int
seed
,
const
double
temp
,
MPI_Comm
root2root
)
:
_lmp
(
lmp
),
inter_comm
(
root2root
)
{
if
(
cvm
::
debug
())
log
(
"Initializing the colvars proxy object.
\n
"
);
_random
=
new
LAMMPS_NS
::
RanPark
(
lmp
,
seed
);
first_timestep
=
true
;
system_force_requested
=
false
;
previous_step
=-
1
;
t_target
=
temp
;
do_exit
=
false
;
restart_every
=
0
;
// User-scripted forces are not available in LAMMPS
force_script_defined
=
false
;
have_scripts
=
false
;
// set input restart name and strip the extension, if present
input_prefix_str
=
std
::
string
(
inp_name
?
inp_name
:
""
);
if
(
input_prefix_str
.
rfind
(
".colvars.state"
)
!=
std
::
string
::
npos
)
input_prefix_str
.
erase
(
input_prefix_str
.
rfind
(
".colvars.state"
),
std
::
string
(
".colvars.state"
).
size
());
// output prefix is always given
output_prefix_str
=
std
::
string
(
out_name
);
// not so for restarts
restart_output_prefix_str
=
std
::
string
(
"rest"
);
// check if it is possible to save output configuration
if
((
!
output_prefix_str
.
size
())
&&
(
!
restart_output_prefix_str
.
size
()))
{
fatal_error
(
"Error: neither the final output state file or "
"the output restart file could be defined, exiting.
\n
"
);
}
// try to extract a restart prefix from a potential restart command.
LAMMPS_NS
::
Output
*
outp
=
_lmp
->
output
;
if
((
outp
->
restart_every_single
>
0
)
&&
(
outp
->
restart1
!=
0
))
{
restart_output_prefix_str
=
std
::
string
(
outp
->
restart1
);
}
else
if
((
outp
->
restart_every_double
>
0
)
&&
(
outp
->
restart2a
!=
0
))
{
restart_output_prefix_str
=
std
::
string
(
outp
->
restart2a
);
}
// trim off unwanted stuff from the restart prefix
if
(
restart_output_prefix_str
.
rfind
(
".*"
)
!=
std
::
string
::
npos
)
restart_output_prefix_str
.
erase
(
restart_output_prefix_str
.
rfind
(
".*"
),
2
);
// initialize multi-replica support, if available
if
(
replica_enabled
())
{
MPI_Comm_rank
(
inter_comm
,
&
inter_me
);
MPI_Comm_size
(
inter_comm
,
&
inter_num
);
}
}
void
colvarproxy_lammps
::
init
(
const
char
*
conf_file
)
{
// create the colvarmodule instance
colvars
=
new
colvarmodule
(
this
);
cvm
::
log
(
"Using LAMMPS interface, version "
+
cvm
::
to_str
(
COLVARPROXY_VERSION
)
+
".
\n
"
);
my_angstrom
=
_lmp
->
force
->
angstrom
;
my_boltzmann
=
_lmp
->
force
->
boltz
;
my_timestep
=
_lmp
->
update
->
dt
*
_lmp
->
force
->
femtosecond
;
// TODO move one or more of these to setup() if needed
colvars
->
read_config_file
(
conf_file
);
colvars
->
setup_input
();
colvars
->
setup_output
();
if
(
_lmp
->
update
->
ntimestep
!=
0
)
{
cvm
::
log
(
"Initializing step number as firstTimestep.
\n
"
);
colvars
->
it
=
colvars
->
it_restart
=
_lmp
->
update
->
ntimestep
;
}
if
(
cvm
::
debug
())
{
log
(
"colvars_atoms = "
+
cvm
::
to_str
(
colvars_atoms
)
+
"
\n
"
);
log
(
"colvars_atoms_ncopies = "
+
cvm
::
to_str
(
colvars_atoms_ncopies
)
+
"
\n
"
);
log
(
"positions = "
+
cvm
::
to_str
(
positions
)
+
"
\n
"
);
log
(
"applied_forces = "
+
cvm
::
to_str
(
applied_forces
)
+
"
\n
"
);
log
(
cvm
::
line_marker
);
log
(
"Info: done initializing the colvars proxy object.
\n
"
);
}
}
colvarproxy_lammps
::~
colvarproxy_lammps
()
{
delete
_random
;
if
(
colvars
!=
NULL
)
{
colvars
->
write_output_files
();
delete
colvars
;
colvars
=
NULL
;
}
}
// re-initialize data where needed
void
colvarproxy_lammps
::
setup
()
{
my_timestep
=
_lmp
->
update
->
dt
*
_lmp
->
force
->
femtosecond
;
colvars
->
setup
();
}
// trigger colvars computation
double
colvarproxy_lammps
::
compute
()
{
if
(
first_timestep
)
{
first_timestep
=
false
;
}
else
{
// Use the time step number inherited from LAMMPS
if
(
_lmp
->
update
->
ntimestep
-
previous_step
==
1
)
colvars
->
it
++
;
// Other cases could mean:
// - run 0
// - beginning of a new run statement
// then the internal counter should not be incremented
}
previous_step
=
_lmp
->
update
->
ntimestep
;
if
(
cvm
::
debug
())
{
cvm
::
log
(
cvm
::
line_marker
+
"colvarproxy_lammps, step no. "
+
cvm
::
to_str
(
colvars
->
it
)
+
"
\n
"
+
"Updating internal data.
\n
"
);
}
// backup applied forces if necessary to calculate system forces
if
(
system_force_requested
)
previous_applied_forces
=
applied_forces
;
// zero the forces on the atoms, so that they can be accumulated by the colvars
for
(
size_t
i
=
0
;
i
<
applied_forces
.
size
();
i
++
)
{
applied_forces
[
i
].
x
=
applied_forces
[
i
].
y
=
applied_forces
[
i
].
z
=
0.0
;
}
bias_energy
=
0.0
;
// call the collective variable module
colvars
->
calc
();
#if 0
for (int i=0; i < colvars_atoms.size(); ++i) {
fprintf(stderr,"CV: atom %d/%d/%d pos: %g %g %g for: %g %g %g\n",
colvars_atoms[i], colvars_atoms_ncopies[i],
positions[i].type, positions[i].x, positions[i].y, positions[i].z,
applied_forces[i].x, applied_forces[i].y, applied_forces[i].z);
}
#endif
return
bias_energy
;
}
void
colvarproxy_lammps
::
serialize_status
(
std
::
string
&
rst
)
{
std
::
ostringstream
os
;
colvars
->
write_restart
(
os
);
rst
=
os
.
str
();
}
// set status from string
bool
colvarproxy_lammps
::
deserialize_status
(
std
::
string
&
rst
)
{
std
::
istringstream
is
;
is
.
str
(
rst
);
if
(
!
colvars
->
read_restart
(
is
))
{
return
false
;
}
else
{
return
true
;
}
}
cvm
::
rvector
colvarproxy_lammps
::
position_distance
(
cvm
::
atom_pos
const
&
pos1
,
cvm
::
atom_pos
const
&
pos2
)
{
double
xtmp
=
pos2
.
x
-
pos1
.
x
;
double
ytmp
=
pos2
.
y
-
pos1
.
y
;
double
ztmp
=
pos2
.
z
-
pos1
.
z
;
_lmp
->
domain
->
minimum_image
(
xtmp
,
ytmp
,
ztmp
);
return
cvm
::
rvector
(
xtmp
,
ytmp
,
ztmp
);
}
cvm
::
real
colvarproxy_lammps
::
position_dist2
(
cvm
::
atom_pos
const
&
pos1
,
cvm
::
atom_pos
const
&
pos2
)
{
double
xtmp
=
pos2
.
x
-
pos1
.
x
;
double
ytmp
=
pos2
.
y
-
pos1
.
y
;
double
ztmp
=
pos2
.
z
-
pos1
.
z
;
_lmp
->
domain
->
minimum_image
(
xtmp
,
ytmp
,
ztmp
);
return
cvm
::
real
(
xtmp
*
xtmp
+
ytmp
*
ytmp
+
ztmp
*
ztmp
);
}
inline
void
colvarproxy_lammps
::
select_closest_image
(
cvm
::
atom_pos
&
pos
,
cvm
::
atom_pos
const
&
ref
)
{
double
xtmp
=
pos
.
x
-
ref
.
x
;
double
ytmp
=
pos
.
y
-
ref
.
y
;
double
ztmp
=
pos
.
z
-
ref
.
z
;
_lmp
->
domain
->
minimum_image
(
xtmp
,
ytmp
,
ztmp
);
pos
.
x
=
ref
.
x
+
xtmp
;
pos
.
y
=
ref
.
y
+
ytmp
;
pos
.
z
=
ref
.
z
+
ztmp
;
}
void
colvarproxy_lammps
::
log
(
std
::
string
const
&
message
)
{
std
::
istringstream
is
(
message
);
std
::
string
line
;
while
(
std
::
getline
(
is
,
line
))
{
if
(
_lmp
->
screen
)
fprintf
(
_lmp
->
screen
,
"colvars: %s
\n
"
,
line
.
c_str
());
if
(
_lmp
->
logfile
)
fprintf
(
_lmp
->
logfile
,
"colvars: %s
\n
"
,
line
.
c_str
());
}
}
void
colvarproxy_lammps
::
error
(
std
::
string
const
&
message
)
{
// In LAMMPS, all errors are fatal
fatal_error
(
message
);
}
void
colvarproxy_lammps
::
fatal_error
(
std
::
string
const
&
message
)
{
log
(
message
);
if
(
!
cvm
::
debug
())
log
(
"If this error message is unclear, try recompiling the "
"colvars library and LAMMPS with -DCOLVARS_DEBUG.
\n
"
);
_lmp
->
error
->
one
(
FLERR
,
"Fatal error in the collective variables module.
\n
"
);
}
void
colvarproxy_lammps
::
exit
(
std
::
string
const
&
message
)
{
log
(
message
);
log
(
"Request to exit the simulation made.
\n
"
);
do_exit
=
true
;
}
enum
e_pdb_field
{
e_pdb_none
,
e_pdb_occ
,
e_pdb_beta
,
e_pdb_x
,
e_pdb_y
,
e_pdb_z
,
e_pdb_ntot
};
e_pdb_field
pdb_field_str2enum
(
std
::
string
const
&
pdb_field_str
)
{
e_pdb_field
pdb_field
=
e_pdb_none
;
if
(
colvarparse
::
to_lower_cppstr
(
pdb_field_str
)
==
colvarparse
::
to_lower_cppstr
(
"O"
))
{
pdb_field
=
e_pdb_occ
;
}
if
(
colvarparse
::
to_lower_cppstr
(
pdb_field_str
)
==
colvarparse
::
to_lower_cppstr
(
"B"
))
{
pdb_field
=
e_pdb_beta
;
}
if
(
colvarparse
::
to_lower_cppstr
(
pdb_field_str
)
==
colvarparse
::
to_lower_cppstr
(
"X"
))
{
pdb_field
=
e_pdb_x
;
}
if
(
colvarparse
::
to_lower_cppstr
(
pdb_field_str
)
==
colvarparse
::
to_lower_cppstr
(
"Y"
))
{
pdb_field
=
e_pdb_y
;
}
if
(
colvarparse
::
to_lower_cppstr
(
pdb_field_str
)
==
colvarparse
::
to_lower_cppstr
(
"Z"
))
{
pdb_field
=
e_pdb_z
;
}
if
(
pdb_field
==
e_pdb_none
)
{
cvm
::
fatal_error
(
"Error: unsupported PDB field,
\"
"
+
pdb_field_str
+
"
\"
.
\n
"
);
}
return
pdb_field
;
}
int
colvarproxy_lammps
::
load_coords
(
char
const
*
pdb_filename
,
std
::
vector
<
cvm
::
atom_pos
>
&
pos
,
const
std
::
vector
<
int
>
&
indices
,
std
::
string
const
&
pdb_field_str
,
double
const
pdb_field_value
)
{
cvm
::
fatal_error
(
"Reading collective variable coordinates "
"from a PDB file is currently not supported.
\n
"
);
return
COLVARS_ERROR
;
}
int
colvarproxy_lammps
::
load_atoms
(
char
const
*
pdb_filename
,
std
::
vector
<
cvm
::
atom
>
&
atoms
,
std
::
string
const
&
pdb_field_str
,
double
const
pdb_field_value
)
{
cvm
::
fatal_error
(
"Selecting collective variable atoms "
"from a PDB file is currently not supported.
\n
"
);
return
COLVARS_ERROR
;
}
int
colvarproxy_lammps
::
backup_file
(
char
const
*
filename
)
{
if
(
std
::
string
(
filename
).
rfind
(
std
::
string
(
".colvars.state"
))
!=
std
::
string
::
npos
)
{
return
my_backup_file
(
filename
,
".old"
);
}
else
{
return
my_backup_file
(
filename
,
".BAK"
);
}
}
int
colvarproxy_lammps
::
init_lammps_atom
(
const
int
&
aid
,
cvm
::
atom
*
atom
)
{
atom
->
id
=
aid
;
atom
->
mass
=
0.0
;
for
(
size_t
i
=
0
;
i
<
colvars_atoms
.
size
();
i
++
)
{
if
(
colvars_atoms
[
i
]
==
aid
)
{
// this atom id was already recorded
colvars_atoms_ncopies
[
i
]
+=
1
;
return
i
;
}
}
// allocate a new slot for this atom
colvars_atoms_ncopies
.
push_back
(
1
);
colvars_atoms
.
push_back
(
aid
);
struct
commdata
c
;
c
.
tag
=
aid
;
c
.
type
=
0
;
c
.
x
=
c
.
y
=
c
.
z
=
0.0
;
positions
.
push_back
(
c
);
total_forces
.
push_back
(
c
);
applied_forces
.
push_back
(
c
);
return
colvars_atoms
.
size
()
-
1
;
}
// multi-replica support
void
colvarproxy_lammps
::
replica_comm_barrier
()
{
MPI_Barrier
(
inter_comm
);
}
int
colvarproxy_lammps
::
replica_comm_recv
(
char
*
msg_data
,
int
buf_len
,
int
src_rep
)
{
MPI_Status
status
;
int
retval
;
retval
=
MPI_Recv
(
msg_data
,
buf_len
,
MPI_CHAR
,
src_rep
,
0
,
inter_comm
,
&
status
);
if
(
retval
==
MPI_SUCCESS
)
{
MPI_Get_count
(
&
status
,
MPI_CHAR
,
&
retval
);
}
else
retval
=
0
;
return
retval
;
}
int
colvarproxy_lammps
::
replica_comm_send
(
char
*
msg_data
,
int
msg_len
,
int
dest_rep
)
{
int
retval
;
retval
=
MPI_Send
(
msg_data
,
msg_len
,
MPI_CHAR
,
dest_rep
,
0
,
inter_comm
);
if
(
retval
==
MPI_SUCCESS
)
{
retval
=
msg_len
;
}
else
retval
=
0
;
return
retval
;
}
// atom member functions, LAMMPS specific implementations
cvm
::
atom
::
atom
(
const
int
&
id
)
{
if
(
cvm
::
debug
())
cvm
::
log
(
"Adding atom "
+
cvm
::
to_str
(
id
)
+
" for collective variables calculation.
\n
"
);
if
(
id
<
0
)
cvm
::
fatal_error
(
"Error: invalid atom ID specified, "
+
cvm
::
to_str
(
id
)
+
"
\n
"
);
int
idx
=
((
colvarproxy_lammps
*
)
cvm
::
proxy
)
->
init_lammps_atom
(
id
,
this
);
if
(
idx
<
0
)
cvm
::
fatal_error
(
"Error: atom ID , "
+
cvm
::
to_str
(
id
)
+
" does not exist.
\n
"
);
this
->
index
=
idx
;
if
(
cvm
::
debug
())
cvm
::
log
(
"The index of this atom in the colvarproxy_lammps arrays is "
+
cvm
::
to_str
(
this
->
index
)
+
".
\n
"
);
this
->
reset_data
();
}
/// For AMBER topologies, the segment id is automatically set to
/// "MAIN" (the segment id assigned by NAMD's AMBER topology parser),
/// and is therefore optional when an AMBER topology is used
cvm
::
atom
::
atom
(
cvm
::
residue_id
const
&
residue
,
std
::
string
const
&
atom_name
,
std
::
string
const
&
segment_id
)
{
cvm
::
fatal_error
(
"Creating collective variable atoms "
"from a PDB file is currently not supported.
\n
"
);
}
// copy constructor
cvm
::
atom
::
atom
(
cvm
::
atom
const
&
a
)
:
index
(
a
.
index
),
id
(
a
.
id
),
mass
(
a
.
mass
)
{
// init_lammps_atom() has already been called by a's constructor, no
// need to call it again
// need to increment the counter anyway
colvarproxy_lammps
*
cp
=
(
colvarproxy_lammps
*
)
cvm
::
proxy
;
cp
->
colvars_atoms_ncopies
[
this
->
index
]
+=
1
;
}
cvm
::
atom
::~
atom
()
{
if
(
this
->
index
>=
0
)
{
colvarproxy_lammps
*
cp
=
(
colvarproxy_lammps
*
)
cvm
::
proxy
;
if
(
cp
->
colvars_atoms_ncopies
[
this
->
index
]
>
0
)
cp
->
colvars_atoms_ncopies
[
this
->
index
]
-=
1
;
}
}
void
cvm
::
atom
::
read_position
()
{
colvarproxy_lammps
const
*
const
cp
=
(
colvarproxy_lammps
*
)
cvm
::
proxy
;
this
->
pos
.
x
=
cp
->
positions
[
this
->
index
].
x
;
this
->
pos
.
y
=
cp
->
positions
[
this
->
index
].
y
;
this
->
pos
.
z
=
cp
->
positions
[
this
->
index
].
z
;
this
->
mass
=
cp
->
positions
[
this
->
index
].
m
;
}
void
cvm
::
atom
::
read_velocity
()
{
cvm
::
fatal_error
(
"Error: read_velocity is not yet implemented.
\n
"
);
}
void
cvm
::
atom
::
read_system_force
()
{
colvarproxy_lammps
const
*
const
cp
=
(
colvarproxy_lammps
*
)
cvm
::
proxy
;
this
->
system_force
.
x
=
cp
->
total_forces
[
this
->
index
].
x
-
cp
->
previous_applied_forces
[
this
->
index
].
x
;
this
->
system_force
.
y
=
cp
->
total_forces
[
this
->
index
].
y
-
cp
->
previous_applied_forces
[
this
->
index
].
y
;
this
->
system_force
.
z
=
cp
->
total_forces
[
this
->
index
].
z
-
cp
->
previous_applied_forces
[
this
->
index
].
z
;
}
void
cvm
::
atom
::
apply_force
(
cvm
::
rvector
const
&
new_force
)
{
colvarproxy_lammps
*
cp
=
(
colvarproxy_lammps
*
)
cvm
::
proxy
;
cp
->
applied_forces
[
this
->
index
].
x
+=
new_force
.
x
;
cp
->
applied_forces
[
this
->
index
].
y
+=
new_force
.
y
;
cp
->
applied_forces
[
this
->
index
].
z
+=
new_force
.
z
;
}
Event Timeline
Log In to Comment