Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79585417
fix_pimd.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, Aug 27, 05:44
Size
22 KB
Mime Type
text/x-c
Expires
Thu, Aug 29, 05:44 (2 d)
Engine
blob
Format
Raw Data
Handle
18323497
Attached To
rLAMMPS lammps
fix_pimd.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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Package FixPIMD
Purpose Quantum Path Integral Algorithm for Quantum Chemistry
Copyright Voth Group @ University of Chicago
Authors Chris Knight & Yuxing Peng (yuxing at uchicago.edu)
Updated Oct-01-2011
Version 1.0
------------------------------------------------------------------------- */
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "fix_pimd.h"
#include "universe.h"
#include "comm.h"
#include "force.h"
#include "atom.h"
#include "domain.h"
#include "update.h"
#include "memory.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
using
namespace
FixConst
;
enum
{
PIMD
,
NMPIMD
,
CMD
};
/* ---------------------------------------------------------------------- */
FixPIMD
::
FixPIMD
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
Fix
(
lmp
,
narg
,
arg
)
{
method
=
PIMD
;
fmass
=
1.0
;
nhc_temp
=
298.15
;
nhc_nchain
=
2
;
sp
=
1.0
;
for
(
int
i
=
3
;
i
<
narg
-
1
;
i
+=
2
)
{
if
(
strcmp
(
arg
[
i
],
"method"
)
==
0
)
{
if
(
strcmp
(
arg
[
i
+
1
],
"pimd"
)
==
0
)
method
=
PIMD
;
else
if
(
strcmp
(
arg
[
i
+
1
],
"nmpimd"
)
==
0
)
method
=
NMPIMD
;
else
if
(
strcmp
(
arg
[
i
+
1
],
"cmd"
)
==
0
)
method
=
CMD
;
else
error
->
universe_all
(
FLERR
,
"Unkown method parameter for fix pimd"
);
}
else
if
(
strcmp
(
arg
[
i
],
"fmass"
)
==
0
)
{
fmass
=
atof
(
arg
[
i
+
1
]);
if
(
fmass
<
0.0
||
fmass
>
1.0
)
error
->
universe_all
(
FLERR
,
"Invalid fmass value for fix pimd"
);
}
else
if
(
strcmp
(
arg
[
i
],
"sp"
)
==
0
)
{
sp
=
atof
(
arg
[
i
+
1
]);
if
(
fmass
<
0.0
)
error
->
universe_all
(
FLERR
,
"Invalid sp value for fix pimd"
);
}
else
if
(
strcmp
(
arg
[
i
],
"temp"
)
==
0
)
{
nhc_temp
=
atof
(
arg
[
i
+
1
]);
if
(
nhc_temp
<
0.0
)
error
->
universe_all
(
FLERR
,
"Invalid temp value for fix pimd"
);
}
else
if
(
strcmp
(
arg
[
i
],
"nhc"
)
==
0
)
{
nhc_nchain
=
atoi
(
arg
[
i
+
1
]);
if
(
nhc_nchain
<
2
)
error
->
universe_all
(
FLERR
,
"Invalid nhc value for fix pimd"
);
}
else
error
->
universe_all
(
arg
[
i
],
i
+
1
,
"Unkown keyword for fix pimd"
);
}
/* Initiation */
max_nsend
=
0
;
tag_send
=
NULL
;
buf_send
=
NULL
;
max_nlocal
=
0
;
buf_recv
=
NULL
;
buf_beads
=
NULL
;
size_plan
=
0
;
plan_send
=
plan_recv
=
NULL
;
M_x2xp
=
M_xp2x
=
M_f2fp
=
M_fp2f
=
NULL
;
lam
=
NULL
;
mode_index
=
NULL
;
mass
=
NULL
;
array_atom
=
NULL
;
nhc_eta
=
NULL
;
nhc_eta_dot
=
NULL
;
nhc_eta_dotdot
=
NULL
;
nhc_eta_mass
=
NULL
;
size_peratom_cols
=
12
*
nhc_nchain
+
3
;
nhc_offset_one_1
=
3
*
nhc_nchain
;
nhc_offset_one_2
=
3
*
nhc_nchain
+
3
;
nhc_size_one_1
=
sizeof
(
double
)
*
nhc_offset_one_1
;
nhc_size_one_2
=
sizeof
(
double
)
*
nhc_offset_one_2
;
restart_peratom
=
1
;
peratom_flag
=
1
;
peratom_freq
=
1
;
global_freq
=
1
;
thermo_energy
=
1
;
vector_flag
=
1
;
size_vector
=
2
;
extvector
=
1
;
comm_forward
=
3
;
atom
->
add_callback
(
0
);
// Call LAMMPS to allocate memory for per-atom array
atom
->
add_callback
(
1
);
// Call LAMMPS to re-assign restart-data for per-atom array
grow_arrays
(
atom
->
nmax
);
// some initilizations
nhc_ready
=
false
;
}
/* ---------------------------------------------------------------------- */
int
FixPIMD
::
setmask
()
{
int
mask
=
0
;
mask
|=
POST_FORCE
;
mask
|=
INITIAL_INTEGRATE
;
mask
|=
FINAL_INTEGRATE
;
return
mask
;
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
init
()
{
if
(
atom
->
map_style
==
0
)
error
->
all
(
FLERR
,
"Fix pimd requires an atom map, see atom_modify"
);
if
(
universe
->
me
==
0
&&
screen
)
fprintf
(
screen
,
"Fix pimd initializing Path-Integral ...
\n
"
);
// prepare the constants
np
=
universe
->
nworlds
;
inverse_np
=
1.0
/
np
;
/* The first solution for the force constant, using SI units
const double Boltzmann = 1.3806488E-23; // SI unit: J/K
const double Plank = 6.6260755E-34; // SI unit: m^2 kg / s
double hbar = Plank / ( 2.0 * M_PI ) * sp;
double beta = 1.0 / ( Boltzmann * input.nh_temp);
// - P / ( beta^2 * hbar^2) SI unit: s^-2
double _fbond = -1.0 / (beta*beta*hbar*hbar) * input.nbeads;
// convert the units: s^-2 -> (kcal/mol) / (g/mol) / (A^2)
fbond = _fbond * 4.184E+26;
*/
/* The current solution, using LAMMPS internal real units */
const
double
Boltzmann
=
force
->
boltz
;
const
double
Plank
=
force
->
hplanck
;
double
hbar
=
Plank
/
(
2.0
*
M_PI
);
double
beta
=
1.0
/
(
Boltzmann
*
nhc_temp
);
double
_fbond
=
1.0
*
np
/
(
beta
*
beta
*
hbar
*
hbar
)
;
omega_np
=
sqrt
(
np
)
/
(
hbar
*
beta
)
*
sqrt
(
force
->
mvv2e
);
fbond
=
-
_fbond
*
force
->
mvv2e
;
if
(
universe
->
me
==
0
)
printf
(
"Fix pimd -P/(beta^2 * hbar^2) = %20.7lE (kcal/mol/A^2)
\n\n
"
,
fbond
);
dtv
=
update
->
dt
;
dtf
=
0.5
*
update
->
dt
*
force
->
ftm2v
;
comm_init
();
mass
=
new
double
[
atom
->
ntypes
+
1
];
if
(
method
==
CMD
||
method
==
NMPIMD
)
nmpimd_init
();
else
for
(
int
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
mass
[
i
]
=
atom
->
mass
[
i
]
/
np
*
fmass
;
if
(
!
nhc_ready
)
nhc_init
();
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
setup
(
int
vflag
)
{
if
(
universe
->
me
==
0
&&
screen
)
fprintf
(
screen
,
"Setting up Path-Integral ...
\n
"
);
post_force
(
vflag
);
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
initial_integrate
(
int
vflag
)
{
nhc_update_v
();
nhc_update_x
();
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
final_integrate
()
{
nhc_update_v
();
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
post_force
(
int
flag
)
{
for
(
int
i
=
0
;
i
<
atom
->
nlocal
;
i
++
)
for
(
int
j
=
0
;
j
<
3
;
j
++
)
atom
->
f
[
i
][
j
]
/=
np
;
comm_exec
(
atom
->
x
);
spring_force
();
if
(
method
==
CMD
||
method
==
NMPIMD
)
{
/* forward comm for the force on ghost atoms */
nmpimd_fill
(
atom
->
f
);
/* inter-partition comm */
comm_exec
(
atom
->
f
);
/* normal-mode transform */
nmpimd_transform
(
buf_beads
,
atom
->
f
,
M_f2fp
[
universe
->
iworld
]);
}
}
/* ----------------------------------------------------------------------
Nose-Hoover Chains
------------------------------------------------------------------------- */
void
FixPIMD
::
nhc_init
()
{
double
tau
=
1.0
/
omega_np
;
double
KT
=
force
->
boltz
*
nhc_temp
;
double
mass0
=
KT
*
tau
*
tau
;
int
max
=
3
*
atom
->
nlocal
;
for
(
int
i
=
0
;
i
<
max
;
i
++
)
{
for
(
int
ichain
=
0
;
ichain
<
nhc_nchain
;
ichain
++
)
{
nhc_eta
[
i
][
ichain
]
=
0.0
;
nhc_eta_dot
[
i
][
ichain
]
=
0.0
;
nhc_eta_dot
[
i
][
ichain
]
=
0.0
;
nhc_eta_dotdot
[
i
][
ichain
]
=
0.0
;
nhc_eta_mass
[
i
][
ichain
]
=
mass0
;
if
((
method
==
CMD
||
method
==
NMPIMD
)
&&
universe
->
iworld
==
0
)
;
else
nhc_eta_mass
[
i
][
ichain
]
*=
fmass
;
}
nhc_eta_dot
[
i
][
nhc_nchain
]
=
0.0
;
for
(
int
ichain
=
1
;
ichain
<
nhc_nchain
;
ichain
++
)
nhc_eta_dotdot
[
i
][
ichain
]
=
(
nhc_eta_mass
[
i
][
ichain
-
1
]
*
nhc_eta_dot
[
i
][
ichain
-
1
]
*
nhc_eta_dot
[
i
][
ichain
-
1
]
*
force
->
mvv2e
-
KT
)
/
nhc_eta_mass
[
i
][
ichain
];
}
// Zero NH acceleration for CMD
if
(
method
==
CMD
&&
universe
->
iworld
==
0
)
for
(
int
i
=
0
;
i
<
max
;
i
++
)
for
(
int
ichain
=
0
;
ichain
<
nhc_nchain
;
ichain
++
)
nhc_eta_dotdot
[
i
][
ichain
]
=
0.0
;
nhc_ready
=
true
;
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
nhc_update_x
()
{
int
n
=
atom
->
nlocal
;
double
**
x
=
atom
->
x
;
double
**
v
=
atom
->
v
;
if
(
method
==
CMD
||
method
==
NMPIMD
)
{
nmpimd_fill
(
atom
->
v
);
comm_exec
(
atom
->
v
);
/* borrow the space of atom->f to store v in cartisian */
v
=
atom
->
f
;
nmpimd_transform
(
buf_beads
,
v
,
M_xp2x
[
universe
->
iworld
]);
}
for
(
int
i
=
0
;
i
<
n
;
i
++
)
{
x
[
i
][
0
]
+=
dtv
*
v
[
i
][
0
];
x
[
i
][
1
]
+=
dtv
*
v
[
i
][
1
];
x
[
i
][
2
]
+=
dtv
*
v
[
i
][
2
];
}
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
nhc_update_v
()
{
int
n
=
atom
->
nlocal
;
int
*
type
=
atom
->
type
;
double
**
v
=
atom
->
v
;
double
**
f
=
atom
->
f
;
for
(
int
i
=
0
;
i
<
n
;
i
++
)
{
double
dtfm
=
dtf
/
mass
[
type
[
i
]];
v
[
i
][
0
]
+=
dtfm
*
f
[
i
][
0
];
v
[
i
][
1
]
+=
dtfm
*
f
[
i
][
1
];
v
[
i
][
2
]
+=
dtfm
*
f
[
i
][
2
];
}
t_sys
=
0.0
;
if
(
method
==
CMD
&&
universe
->
iworld
==
0
)
return
;
double
expfac
;
int
nmax
=
3
*
atom
->
nlocal
;
double
KT
=
force
->
boltz
*
nhc_temp
;
double
kecurrent
,
t_current
;
double
dthalf
=
0.5
*
update
->
dt
;
double
dt4
=
0.25
*
update
->
dt
;
double
dt8
=
0.125
*
update
->
dt
;
for
(
int
i
=
0
;
i
<
nmax
;
i
++
)
{
int
iatm
=
i
/
3
;
int
idim
=
i
%
3
;
double
*
vv
=
v
[
iatm
];
kecurrent
=
mass
[
type
[
iatm
]]
*
vv
[
idim
]
*
vv
[
idim
]
*
force
->
mvv2e
;
t_current
=
kecurrent
/
force
->
boltz
;
double
*
eta
=
nhc_eta
[
i
];
double
*
eta_dot
=
nhc_eta_dot
[
i
];
double
*
eta_dotdot
=
nhc_eta_dotdot
[
i
];
eta_dotdot
[
0
]
=
(
kecurrent
-
KT
)
/
nhc_eta_mass
[
i
][
0
];
for
(
int
ichain
=
nhc_nchain
-
1
;
ichain
>
0
;
ichain
--
)
{
expfac
=
exp
(
-
dt8
*
eta_dot
[
ichain
+
1
]);
eta_dot
[
ichain
]
*=
expfac
;
eta_dot
[
ichain
]
+=
eta_dotdot
[
ichain
]
*
dt4
;
eta_dot
[
ichain
]
*=
expfac
;
}
expfac
=
exp
(
-
dt8
*
eta_dot
[
1
]);
eta_dot
[
0
]
*=
expfac
;
eta_dot
[
0
]
+=
eta_dotdot
[
0
]
*
dt4
;
eta_dot
[
0
]
*=
expfac
;
// Update particle velocities half-step
double
factor_eta
=
exp
(
-
dthalf
*
eta_dot
[
0
]);
vv
[
idim
]
*=
factor_eta
;
t_current
*=
(
factor_eta
*
factor_eta
);
kecurrent
=
force
->
boltz
*
t_current
;
eta_dotdot
[
0
]
=
(
kecurrent
-
KT
)
/
nhc_eta_mass
[
i
][
0
];
for
(
int
ichain
=
0
;
ichain
<
nhc_nchain
;
ichain
++
)
eta
[
ichain
]
+=
dthalf
*
eta_dot
[
ichain
];
eta_dot
[
0
]
*=
expfac
;
eta_dot
[
0
]
+=
eta_dotdot
[
0
]
*
dt4
;
eta_dot
[
0
]
*=
expfac
;
for
(
int
ichain
=
1
;
ichain
<
nhc_nchain
;
ichain
++
)
{
expfac
=
exp
(
-
dt8
*
eta_dot
[
ichain
+
1
]);
eta_dot
[
ichain
]
*=
expfac
;
eta_dotdot
[
ichain
]
=
(
nhc_eta_mass
[
i
][
ichain
-
1
]
*
eta_dot
[
ichain
-
1
]
*
eta_dot
[
ichain
-
1
]
-
KT
)
/
nhc_eta_mass
[
i
][
ichain
];
eta_dot
[
ichain
]
+=
eta_dotdot
[
ichain
]
*
dt4
;
eta_dot
[
ichain
]
*=
expfac
;
}
t_sys
+=
t_current
;
}
t_sys
/=
nmax
;
}
/* ----------------------------------------------------------------------
Normal Mode PIMD
------------------------------------------------------------------------- */
void
FixPIMD
::
nmpimd_init
()
{
memory
->
create
(
M_x2xp
,
np
,
np
,
"fix_feynman:M_x2xp"
);
memory
->
create
(
M_xp2x
,
np
,
np
,
"fix_feynman:M_xp2x"
);
memory
->
create
(
M_f2fp
,
np
,
np
,
"fix_feynman:M_f2fp"
);
memory
->
create
(
M_fp2f
,
np
,
np
,
"fix_feynman:M_fp2f"
);
lam
=
(
double
*
)
memory
->
smalloc
(
sizeof
(
double
)
*
np
,
"FixPIMD::lam"
);
// Set up eigenvalues
lam
[
0
]
=
0.0
;
if
(
np
%
2
==
0
)
lam
[
np
-
1
]
=
4.0
*
np
;
for
(
int
i
=
2
;
i
<=
np
/
2
;
i
++
)
{
lam
[
2
*
i
-
3
]
=
lam
[
2
*
i
-
2
]
=
2.0
*
np
*
(
1.0
-
1.0
*
cos
(
2.0
*
M_PI
*
(
i
-
1
)
/
np
));
}
// Set up eigenvectors for non-degenerated modes
for
(
int
i
=
0
;
i
<
np
;
i
++
)
{
M_x2xp
[
0
][
i
]
=
1.0
/
np
;
if
(
np
%
2
==
0
)
M_x2xp
[
np
-
1
][
i
]
=
1.0
/
np
*
pow
(
-
1.0
,
i
);
}
// Set up eigenvectors for degenerated modes
for
(
int
i
=
0
;
i
<
(
np
-
1
)
/
2
;
i
++
)
for
(
int
j
=
0
;
j
<
np
;
j
++
)
{
M_x2xp
[
2
*
i
+
1
][
j
]
=
sqrt
(
2.0
)
*
cos
(
2.0
*
M_PI
*
(
i
+
1
)
*
j
/
np
)
/
np
;
M_x2xp
[
2
*
i
+
2
][
j
]
=
-
sqrt
(
2.0
)
*
sin
(
2.0
*
M_PI
*
(
i
+
1
)
*
j
/
np
)
/
np
;
}
// Set up Ut
for
(
int
i
=
0
;
i
<
np
;
i
++
)
for
(
int
j
=
0
;
j
<
np
;
j
++
)
{
M_xp2x
[
i
][
j
]
=
M_x2xp
[
j
][
i
]
*
np
;
M_f2fp
[
i
][
j
]
=
M_x2xp
[
i
][
j
]
*
np
;
M_fp2f
[
i
][
j
]
=
M_xp2x
[
i
][
j
];
}
// Set up masses
int
iworld
=
universe
->
iworld
;
for
(
int
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
{
mass
[
i
]
=
atom
->
mass
[
i
];
if
(
iworld
)
{
mass
[
i
]
*=
lam
[
iworld
];
mass
[
i
]
*=
fmass
;
}
}
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
nmpimd_fill
(
double
**
ptr
)
{
comm_ptr
=
ptr
;
comm
->
forward_comm_fix
(
this
);
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
nmpimd_transform
(
double
**
src
,
double
**
des
,
double
*
vector
)
{
int
n
=
atom
->
nlocal
;
int
m
=
0
;
for
(
int
i
=
0
;
i
<
n
;
i
++
)
for
(
int
d
=
0
;
d
<
3
;
d
++
)
{
des
[
i
][
d
]
=
0.0
;
for
(
int
j
=
0
;
j
<
np
;
j
++
)
{
des
[
i
][
d
]
+=
(
src
[
j
][
m
]
*
vector
[
j
]);
}
m
++
;
}
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
spring_force
()
{
spring_energy
=
0.0
;
double
**
x
=
atom
->
x
;
double
**
f
=
atom
->
f
;
double
*
_mass
=
atom
->
mass
;
int
*
type
=
atom
->
type
;
int
nlocal
=
atom
->
nlocal
;
double
*
xlast
=
buf_beads
[
x_last
];
double
*
xnext
=
buf_beads
[
x_next
];
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
double
delx1
=
xlast
[
0
]
-
x
[
i
][
0
];
double
dely1
=
xlast
[
1
]
-
x
[
i
][
1
];
double
delz1
=
xlast
[
2
]
-
x
[
i
][
2
];
xlast
+=
3
;
domain
->
minimum_image
(
delx1
,
dely1
,
delz1
);
double
delx2
=
xnext
[
0
]
-
x
[
i
][
0
];
double
dely2
=
xnext
[
1
]
-
x
[
i
][
1
];
double
delz2
=
xnext
[
2
]
-
x
[
i
][
2
];
xnext
+=
3
;
domain
->
minimum_image
(
delx2
,
dely2
,
delz2
);
double
ff
=
fbond
*
_mass
[
type
[
i
]];
double
dx
=
delx1
+
delx2
;
double
dy
=
dely1
+
dely2
;
double
dz
=
delz1
+
delz2
;
f
[
i
][
0
]
-=
(
dx
)
*
ff
;
f
[
i
][
1
]
-=
(
dy
)
*
ff
;
f
[
i
][
2
]
-=
(
dz
)
*
ff
;
spring_energy
+=
(
dx
*
dx
+
dy
*
dy
+
dz
*
dz
);
}
}
/* ----------------------------------------------------------------------
Comm operations
------------------------------------------------------------------------- */
void
FixPIMD
::
comm_init
()
{
if
(
size_plan
)
{
delete
[]
plan_send
;
delete
[]
plan_recv
;
}
if
(
method
==
PIMD
)
{
size_plan
=
2
;
plan_send
=
new
int
[
2
];
plan_recv
=
new
int
[
2
];
mode_index
=
new
int
[
2
];
int
rank_last
=
universe
->
me
-
comm
->
nprocs
;
int
rank_next
=
universe
->
me
+
comm
->
nprocs
;
if
(
rank_last
<
0
)
rank_last
+=
universe
->
nprocs
;
if
(
rank_next
>=
universe
->
nprocs
)
rank_next
-=
universe
->
nprocs
;
plan_send
[
0
]
=
rank_next
;
plan_send
[
1
]
=
rank_last
;
plan_recv
[
0
]
=
rank_last
;
plan_recv
[
1
]
=
rank_next
;
mode_index
[
0
]
=
0
;
mode_index
[
1
]
=
1
;
x_last
=
1
;
x_next
=
0
;
}
else
{
size_plan
=
np
-
1
;
plan_send
=
new
int
[
size_plan
];
plan_recv
=
new
int
[
size_plan
];
mode_index
=
new
int
[
size_plan
];
for
(
int
i
=
0
;
i
<
size_plan
;
i
++
)
{
plan_send
[
i
]
=
universe
->
me
+
comm
->
nprocs
*
(
i
+
1
);
if
(
plan_send
[
i
]
>=
universe
->
nprocs
)
plan_send
[
i
]
-=
universe
->
nprocs
;
plan_recv
[
i
]
=
universe
->
me
-
comm
->
nprocs
*
(
i
+
1
);
if
(
plan_recv
[
i
]
<
0
)
plan_recv
[
i
]
+=
universe
->
nprocs
;
mode_index
[
i
]
=
(
universe
->
iworld
+
i
+
1
)
%
(
universe
->
nworlds
);
}
x_next
=
(
universe
->
iworld
+
1
+
universe
->
nworlds
)
%
(
universe
->
nworlds
);
x_last
=
(
universe
->
iworld
-
1
+
universe
->
nworlds
)
%
(
universe
->
nworlds
);
}
if
(
buf_beads
)
{
for
(
int
i
=
0
;
i
<
np
;
i
++
)
if
(
buf_beads
[
i
])
delete
[]
buf_beads
[
i
];
delete
[]
buf_beads
;
}
buf_beads
=
new
double
*
[
np
];
for
(
int
i
=
0
;
i
<
np
;
i
++
)
buf_beads
[
i
]
=
NULL
;
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
comm_exec
(
double
**
ptr
)
{
int
nlocal
=
atom
->
nlocal
;
if
(
nlocal
>
max_nlocal
)
{
max_nlocal
=
nlocal
+
200
;
int
size
=
sizeof
(
double
)
*
max_nlocal
*
3
;
buf_recv
=
(
double
*
)
memory
->
srealloc
(
buf_recv
,
size
,
"FixPIMD:x_recv"
);
for
(
int
i
=
0
;
i
<
np
;
i
++
)
buf_beads
[
i
]
=
(
double
*
)
memory
->
srealloc
(
buf_beads
[
i
],
size
,
"FixPIMD:x_beads[i]"
);
}
// copy local positions
memcpy
(
buf_beads
[
universe
->
iworld
],
&
(
ptr
[
0
][
0
]),
sizeof
(
double
)
*
nlocal
*
3
);
// go over comm plans
for
(
int
iplan
=
0
;
iplan
<
size_plan
;
iplan
++
)
{
// sendrecv nlocal
int
nsend
;
MPI_Sendrecv
(
&
(
nlocal
),
1
,
MPI_INT
,
plan_send
[
iplan
],
0
,
&
(
nsend
),
1
,
MPI_INT
,
plan_recv
[
iplan
],
0
,
universe
->
uworld
,
MPI_STATUS_IGNORE
);
// allocate arrays
if
(
nsend
>
max_nsend
)
{
max_nsend
=
nsend
+
200
;
tag_send
=
(
int
*
)
memory
->
srealloc
(
tag_send
,
sizeof
(
int
)
*
max_nsend
,
"FixPIMD:tag_send"
);
buf_send
=
(
double
*
)
memory
->
srealloc
(
buf_send
,
sizeof
(
double
)
*
max_nsend
*
3
,
"FixPIMD:x_send"
);
}
// send tags
MPI_Sendrecv
(
atom
->
tag
,
nlocal
,
MPI_INT
,
plan_send
[
iplan
],
0
,
tag_send
,
nsend
,
MPI_INT
,
plan_recv
[
iplan
],
0
,
universe
->
uworld
,
MPI_STATUS_IGNORE
);
// wrap positions
double
*
wrap_ptr
=
buf_send
;
int
ncpy
=
sizeof
(
double
)
*
3
;
for
(
int
i
=
0
;
i
<
nsend
;
i
++
)
{
int
index
=
atom
->
map
(
tag_send
[
i
]);
if
(
index
<
0
)
{
char
error_line
[
256
];
sprintf
(
error_line
,
"Atom "
TAGINT_FORMAT
" is missing at world [%d] "
"rank [%d] required by rank [%d] ("
TAGINT_FORMAT
", "
TAGINT_FORMAT
", "
TAGINT_FORMAT
").
\n
"
,
tag_send
[
i
],
universe
->
iworld
,
comm
->
me
,
plan_recv
[
iplan
],
atom
->
tag
[
0
],
atom
->
tag
[
1
],
atom
->
tag
[
2
]);
error
->
universe_one
(
FLERR
,
error_line
);
}
memcpy
(
wrap_ptr
,
ptr
[
index
],
ncpy
);
wrap_ptr
+=
3
;
}
// sendrecv x
MPI_Sendrecv
(
buf_send
,
nsend
*
3
,
MPI_DOUBLE
,
plan_recv
[
iplan
],
0
,
buf_recv
,
nlocal
*
3
,
MPI_DOUBLE
,
plan_send
[
iplan
],
0
,
universe
->
uworld
,
MPI_STATUS_IGNORE
);
// copy x
memcpy
(
buf_beads
[
mode_index
[
iplan
]],
buf_recv
,
sizeof
(
double
)
*
nlocal
*
3
);
}
}
/* ---------------------------------------------------------------------- */
int
FixPIMD
::
pack_forward_comm
(
int
n
,
int
*
list
,
double
*
buf
,
int
pbc_flag
,
int
*
pbc
)
{
int
i
,
j
,
m
;
m
=
0
;
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
buf
[
m
++
]
=
comm_ptr
[
j
][
0
];
buf
[
m
++
]
=
comm_ptr
[
j
][
1
];
buf
[
m
++
]
=
comm_ptr
[
j
][
2
];
}
return
m
;
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
unpack_forward_comm
(
int
n
,
int
first
,
double
*
buf
)
{
int
i
,
m
,
last
;
m
=
0
;
last
=
first
+
n
;
for
(
i
=
first
;
i
<
last
;
i
++
)
{
comm_ptr
[
i
][
0
]
=
buf
[
m
++
];
comm_ptr
[
i
][
1
]
=
buf
[
m
++
];
comm_ptr
[
i
][
2
]
=
buf
[
m
++
];
}
}
/* ----------------------------------------------------------------------
Memory operations
------------------------------------------------------------------------- */
double
FixPIMD
::
memory_usage
()
{
double
bytes
=
0
;
bytes
=
atom
->
nmax
*
size_peratom_cols
*
sizeof
(
double
);
return
bytes
;
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
grow_arrays
(
int
nmax
)
{
if
(
nmax
==
0
)
return
;
int
count
=
nmax
*
3
;
memory
->
grow
(
array_atom
,
nmax
,
size_peratom_cols
,
"FixPIMD::array_atom"
);
memory
->
grow
(
nhc_eta
,
count
,
nhc_nchain
,
"FixPIMD::nh_eta"
);
memory
->
grow
(
nhc_eta_dot
,
count
,
nhc_nchain
+
1
,
"FixPIMD::nh_eta_dot"
);
memory
->
grow
(
nhc_eta_dotdot
,
count
,
nhc_nchain
,
"FixPIMD::nh_eta_dotdot"
);
memory
->
grow
(
nhc_eta_mass
,
count
,
nhc_nchain
,
"FixPIMD::nh_eta_mass"
);
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
copy_arrays
(
int
i
,
int
j
,
int
delflag
)
{
int
i_pos
=
i
*
3
;
int
j_pos
=
j
*
3
;
memcpy
(
nhc_eta
[
j_pos
],
nhc_eta
[
i_pos
],
nhc_size_one_1
);
memcpy
(
nhc_eta_dot
[
j_pos
],
nhc_eta_dot
[
i_pos
],
nhc_size_one_2
);
memcpy
(
nhc_eta_dotdot
[
j_pos
],
nhc_eta_dotdot
[
i_pos
],
nhc_size_one_1
);
memcpy
(
nhc_eta_mass
[
j_pos
],
nhc_eta_mass
[
i_pos
],
nhc_size_one_1
);
}
/* ---------------------------------------------------------------------- */
int
FixPIMD
::
pack_exchange
(
int
i
,
double
*
buf
)
{
int
offset
=
0
;
int
pos
=
i
*
3
;
memcpy
(
buf
+
offset
,
nhc_eta
[
pos
],
nhc_size_one_1
);
offset
+=
nhc_offset_one_1
;
memcpy
(
buf
+
offset
,
nhc_eta_dot
[
pos
],
nhc_size_one_2
);
offset
+=
nhc_offset_one_2
;
memcpy
(
buf
+
offset
,
nhc_eta_dotdot
[
pos
],
nhc_size_one_1
);
offset
+=
nhc_offset_one_1
;
memcpy
(
buf
+
offset
,
nhc_eta_mass
[
pos
],
nhc_size_one_1
);
offset
+=
nhc_offset_one_1
;
return
size_peratom_cols
;
}
/* ---------------------------------------------------------------------- */
int
FixPIMD
::
unpack_exchange
(
int
nlocal
,
double
*
buf
)
{
int
offset
=
0
;
int
pos
=
nlocal
*
3
;
memcpy
(
nhc_eta
[
pos
],
buf
+
offset
,
nhc_size_one_1
);
offset
+=
nhc_offset_one_1
;
memcpy
(
nhc_eta_dot
[
pos
],
buf
+
offset
,
nhc_size_one_2
);
offset
+=
nhc_offset_one_2
;
memcpy
(
nhc_eta_dotdot
[
pos
],
buf
+
offset
,
nhc_size_one_1
);
offset
+=
nhc_offset_one_1
;
memcpy
(
nhc_eta_mass
[
pos
],
buf
+
offset
,
nhc_size_one_1
);
offset
+=
nhc_offset_one_1
;
return
size_peratom_cols
;
}
/* ---------------------------------------------------------------------- */
int
FixPIMD
::
pack_restart
(
int
i
,
double
*
buf
)
{
int
offset
=
0
;
int
pos
=
i
*
3
;
buf
[
offset
++
]
=
size_peratom_cols
+
1
;
memcpy
(
buf
+
offset
,
nhc_eta
[
pos
],
nhc_size_one_1
);
offset
+=
nhc_offset_one_1
;
memcpy
(
buf
+
offset
,
nhc_eta_dot
[
pos
],
nhc_size_one_2
);
offset
+=
nhc_offset_one_2
;
memcpy
(
buf
+
offset
,
nhc_eta_dotdot
[
pos
],
nhc_size_one_1
);
offset
+=
nhc_offset_one_1
;
memcpy
(
buf
+
offset
,
nhc_eta_mass
[
pos
],
nhc_size_one_1
);
offset
+=
nhc_offset_one_1
;
return
size_peratom_cols
+
1
;
}
/* ---------------------------------------------------------------------- */
void
FixPIMD
::
unpack_restart
(
int
nlocal
,
int
nth
)
{
double
**
extra
=
atom
->
extra
;
// skip to Nth set of extra values
int
m
=
0
;
for
(
int
i
=
0
;
i
<
nth
;
i
++
)
m
+=
static_cast
<
int
>
(
extra
[
nlocal
][
m
]);
m
++
;
int
pos
=
nlocal
*
3
;
memcpy
(
nhc_eta
[
pos
],
extra
[
nlocal
]
+
m
,
nhc_size_one_1
);
m
+=
nhc_offset_one_1
;
memcpy
(
nhc_eta_dot
[
pos
],
extra
[
nlocal
]
+
m
,
nhc_size_one_2
);
m
+=
nhc_offset_one_2
;
memcpy
(
nhc_eta_dotdot
[
pos
],
extra
[
nlocal
]
+
m
,
nhc_size_one_1
);
m
+=
nhc_offset_one_1
;
memcpy
(
nhc_eta_mass
[
pos
],
extra
[
nlocal
]
+
m
,
nhc_size_one_1
);
m
+=
nhc_offset_one_1
;
nhc_ready
=
true
;
}
/* ---------------------------------------------------------------------- */
int
FixPIMD
::
maxsize_restart
()
{
return
size_peratom_cols
+
1
;
}
/* ---------------------------------------------------------------------- */
int
FixPIMD
::
size_restart
(
int
nlocal
)
{
return
size_peratom_cols
+
1
;
}
/* ---------------------------------------------------------------------- */
double
FixPIMD
::
compute_vector
(
int
n
)
{
if
(
n
==
0
)
{
return
spring_energy
;
}
if
(
n
==
1
)
{
return
t_sys
;
}
return
0.0
;
}
Event Timeline
Log In to Comment