forked from lijiext/lammps
Axel Kohlmeyer
GitHub
commit
e7dcb79ac5
|
|
@ -150,6 +150,7 @@ OPT.
|
|||
* :doc:`orient/bcc <fix_orient>`
|
||||
* :doc:`orient/fcc <fix_orient>`
|
||||
* :doc:`orient/eco <fix_orient_eco>`
|
||||
* :doc:`pafi <fix_pafi>`
|
||||
* :doc:`phonon <fix_phonon>`
|
||||
* :doc:`pimd <fix_pimd>`
|
||||
* :doc:`planeforce <fix_planeforce>`
|
||||
|
|
|
|||
|
|
@ -1903,6 +1903,12 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
|
|||
*Compute %s does not allow use of dynamic group*
|
||||
Dynamic groups have not yet been enabled for this compute.
|
||||
|
||||
*Compute for fix pafi does not calculate a local array*
|
||||
Self-explanatory.
|
||||
|
||||
*Compute for fix pafi must have 9 fields per atom*
|
||||
Self-explanatory.
|
||||
|
||||
*Compute ID for compute chunk /atom does not exist*
|
||||
Self-explanatory.
|
||||
|
||||
|
|
|
|||
|
|
@ -293,6 +293,7 @@ accelerated styles exist.
|
|||
* :doc:`orient/bcc <fix_orient>` - add grain boundary migration force for BCC
|
||||
* :doc:`orient/fcc <fix_orient>` - add grain boundary migration force for FCC
|
||||
* :doc:`orient/eco <fix_orient_eco>` - add generalized grain boundary migration force
|
||||
* :doc:`pafi <fix_pafi>` - constrained force averages on hyper-planes to compute free energies (PAFI)
|
||||
* :doc:`phonon <fix_phonon>` - calculate dynamical matrix from MD simulations
|
||||
* :doc:`pimd <fix_pimd>` - Feynman path integral molecular dynamics
|
||||
* :doc:`planeforce <fix_planeforce>` - constrain atoms to move in a plane
|
||||
|
|
|
|||
|
|
@ -0,0 +1,99 @@
|
|||
.. index:: fix pafi
|
||||
|
||||
fix pafi command
|
||||
================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID pafi compute-ID Temp Tdamp seed keyword values...
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* pafi = style name of this fix command
|
||||
* compute-ID = ID of a :doc:`compute property/atom <compute_property_atom>` that holds data used by this fix
|
||||
* Temp = desired temperature (temperature units)
|
||||
* Tdamp = damping parameter (time units)
|
||||
* seed = random number seed to use for white noise (positive integer)
|
||||
* keyword = *overdamped* or *com*
|
||||
|
||||
.. parsed-literal::
|
||||
*overdamped* value = *yes* or *no* or 1 or 0
|
||||
*yes* or 1 = Brownian (overdamped) integration in hyperplane
|
||||
*no* or 0 = Langevin integration in hyperplane
|
||||
*com* value = *yes* or *no* or 1 or 0
|
||||
*yes* or 1 = zero linear momentum, fixing center or mass (recommended)
|
||||
*no* or 0 = do not zero linear momentum, allowing center of mass drift
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute pa all property/atom d_nx d_ny d_nz d_dnx d_dny d_dnz d_ddnx d_ddny d_ddnz
|
||||
run 0 post no
|
||||
fix hp all pafi pa 500.0 0.01 434 overdamped yes
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Perform Brownian or Langevin integration whilst constraining the system to lie
|
||||
in some hyperplane, which is expected to be the tangent plane to some reference
|
||||
pathway in a solid state system. The instantaneous value of a modified force
|
||||
projection is also calculated, whose time integral can be shown to be equal to
|
||||
the true free energy gradient along the minimum free energy path local to the reference pathway.
|
||||
A detailed discussion of the projection technique can be found in :ref:`(Swinburne) <Swinburne>`.
|
||||
|
||||
This fix can be used with LAMMPS as demonstrated in examples/USER/misc/pafi,
|
||||
though it is primarily intended to be coupled with the PAFI C++ code, developed
|
||||
at `https://github.com/tomswinburne/pafi <https://github.com/tomswinburne/pafi>`_,
|
||||
which distributes multiple LAMMPS workers in parallel to compute and collate
|
||||
hyperplane-constrained averages, allowing the calculation of free energy barriers and pathways.
|
||||
|
||||
A :doc:`compute property/atom <compute_property_atom>` must be provided with 9 fields per atom coordinate,
|
||||
which in order are the x,y,z coordinates of a configuration on the reference path,
|
||||
the x,y,z coordinates of the path tangent (derivative of path position with path coordinate)
|
||||
and the x,y,z coordinates of the change in tangent (derivative of path tangent with path coordinate).
|
||||
|
||||
A 4-element vector is also calculated by this fix. The 4 components are the
|
||||
modified projected force, its square, the expected projection of the minimum
|
||||
free energy path tangent on the reference path tangent and the minimum image
|
||||
distance between the current configuration and the reference configuration,
|
||||
projected along the path tangent. This latter value should be essentially zero.
|
||||
|
||||
.. note::
|
||||
When com=yes/1, which is recommended, the provided tangent vector must also
|
||||
have zero center of mass. This can be achieved by subtracting from each
|
||||
coordinate of the path tangent the average x,y,z value. The PAFI C++ code
|
||||
(see above) can generate these paths for use in LAMMPS.
|
||||
|
||||
.. note::
|
||||
When overdamped=yes/1, the Tdamp parameter should be around 5-10 times smaller
|
||||
than that used in typical Langevin integration.
|
||||
See :doc:`fix langevin <fix_langevin>` for typical values.
|
||||
|
||||
|
||||
Restart, fix_modify, output, run start/stop, minimize info
|
||||
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
No information about this fix is written to :doc:`binary restart files <restart>`.
|
||||
None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
|
||||
|
||||
This fix produces a global vector each timestep which can be accessed by various :doc:`output commands <Howto_output>`.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This fix is part of the USER-MISC package. It is only enabled if
|
||||
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
|
||||
|
||||
|
||||
Default
|
||||
"""""""
|
||||
The option defaults are com = *yes*, overdamped = *no*
|
||||
|
||||
----------
|
||||
|
||||
.. _Swinburne:
|
||||
|
||||
**(Swinburne)** Swinburne and Marinica, Physical Review Letters, 120, 1 (2018)
|
||||
|
|
@ -1241,6 +1241,7 @@ hydrostatically
|
|||
Hynninen
|
||||
Hyoungki
|
||||
hyperdynamics
|
||||
hyperplane
|
||||
hyperradius
|
||||
hyperspherical
|
||||
hysteretic
|
||||
|
|
@ -1724,6 +1725,7 @@ Maras
|
|||
Marchetti
|
||||
Marchi
|
||||
Mariella
|
||||
Marinica
|
||||
Marrink
|
||||
Marroquin
|
||||
Marsaglia
|
||||
|
|
@ -2272,6 +2274,7 @@ oxRNA
|
|||
packings
|
||||
padua
|
||||
Padua
|
||||
pafi
|
||||
palegoldenrod
|
||||
palegreen
|
||||
paleturquoise
|
||||
|
|
@ -2947,6 +2950,7 @@ sw
|
|||
Swegat
|
||||
swiggle
|
||||
Swiler
|
||||
Swinburne
|
||||
Swol
|
||||
Swope
|
||||
Sx
|
||||
|
|
|
|||
|
|
@ -0,0 +1,2 @@
|
|||
*.dat
|
||||
lmp*
|
||||
|
|
@ -0,0 +1,40 @@
|
|||
Run an example constrained sample for USER-PAFI calculation
|
||||
|
||||
Author:
|
||||
Thomas Swinburne, CNRS / CINaM, Marseille
|
||||
swinburne@cinam.univ-mrs.fr
|
||||
tomswinburne.github.io
|
||||
|
||||
|
||||
PAFI performs constrained force averages on hyperplanes to a reference pathway in order to compute unbiased free energy barriers
|
||||
|
||||
|
||||
Paper:
|
||||
T. D. Swinburne and M. C. Marinica
|
||||
Unsupervised calculation of free energy barriers in large crystalline systems
|
||||
Physical Review Letters 276, 1 2018
|
||||
|
||||
|
||||
Also see https://github.com/tomswinburne/pafi/ for the PAFI code which distributes
|
||||
multiple LAMMPS workers in parallel to compute and collate fix pafi averages,
|
||||
allowing the calculation of free energy barriers and minimum free energy paths.
|
||||
|
||||
To compile:
|
||||
|
||||
make yes-user-misc # for PAFI
|
||||
make yes-manybody # for EAM potential
|
||||
make machine # for binary (machine= e.g. mpi)
|
||||
|
||||
To run the example from this folder:
|
||||
```
|
||||
mpirun -np NPROCS /path/to/lammps/src/lmp_machine -in in.pafi
|
||||
```
|
||||
|
||||
This example executes a hyperplane constrained overdamped Langevin simulation
|
||||
whilst recording the projected free energy gradient, then minimizes in-plane,
|
||||
then removes the hyperplane constraint and minimizes the whole system.
|
||||
If the temperature is too high the system will leave the local minima of
|
||||
the constrained system; PAFI checks for this also.
|
||||
|
||||
|
||||
To generate your own path, please see the LAMMPS documentation for more details
|
||||
|
|
@ -0,0 +1,86 @@
|
|||
units metal
|
||||
atom_style atomic
|
||||
atom_modify map array sort 0 0.0
|
||||
neigh_modify every 2 delay 10 check yes page 1000000 one 100000
|
||||
|
||||
## read in path data using fix property/atom- here 4th image of a NEB
|
||||
fix pa all property/atom d_nx d_ny d_nz d_dnx d_dny d_dnz d_ddnx d_ddny d_ddnz
|
||||
read_data pafipath.4.data fix pa NULL PafiPath
|
||||
|
||||
## EAM potential
|
||||
pair_style eam/fs
|
||||
pair_coeff * * ../../../../potentials/Fe_mm.eam.fs Fe
|
||||
|
||||
mass * 55.85
|
||||
thermo 100
|
||||
min_style fire
|
||||
|
||||
|
||||
compute pe all pe
|
||||
variable pe equal pe
|
||||
run 0
|
||||
print "energy=${pe}"
|
||||
|
||||
## compute property/atom to access relevant fields
|
||||
compute pa all property/atom d_nx d_ny d_nz d_dnx d_dny d_dnz d_ddnx d_ddny d_ddnz
|
||||
run 0
|
||||
|
||||
## fix name group-id pafi compute-id temperature tdamp seed overdamped 0/1 com 0/1
|
||||
fix hp all pafi pa 500.0 0.01 434 overdamped no com yes
|
||||
|
||||
run 0
|
||||
minimize 0 0 1000 1000 # best if using NEB path
|
||||
compute dx all displace/atom
|
||||
compute dmag all reduce max c_dx[4]
|
||||
variable maxjump equal sqrt(c_dmag)
|
||||
|
||||
thermo_style custom step pe temp v_maxjump
|
||||
|
||||
|
||||
variable dn equal f_hp[4]*f_hp[4] # should be zero to float precision
|
||||
|
||||
|
||||
fix af all ave/time 1 1000 1000 f_hp[1] f_hp[2] f_hp[3] v_dn
|
||||
|
||||
variable adn equal sqrt(f_af[4]) # should be very small (approx. 1e-6 A)
|
||||
variable apf equal f_af[1]
|
||||
|
||||
|
||||
run 1000
|
||||
minimize 0 0 1000 1000
|
||||
|
||||
variable s_pe equal ${pe}
|
||||
variable s_apf equal ${apf}
|
||||
variable s_adn equal ${adn}
|
||||
variable s_maxjump equal ${maxjump}
|
||||
|
||||
unfix hp
|
||||
unfix pa
|
||||
unfix af
|
||||
minimize 0 0 1000 1000
|
||||
|
||||
|
||||
print """
|
||||
---------- PAFI RESULTS --------
|
||||
|
||||
During run:
|
||||
Average Distance From Hyperplane = ${s_adn}A (Should be very small, around 1e-5A),
|
||||
Average Force Projection Along MFEP = ${s_apf} eV/A
|
||||
"""
|
||||
print """
|
||||
In-plane minimization post-run:
|
||||
energy = ${s_pe} eV
|
||||
|
||||
"""
|
||||
if "${s_maxjump} > 0.1" then &
|
||||
"print 'Max Atomic Displacement = ${s_maxjump}A > 0.1A'" &
|
||||
"print ' => Possible shallow in-plane metastability: Reduce time in basin and/or decrease tdamp'" &
|
||||
else &
|
||||
"print 'Max Atomic Displacement = ${s_maxjump}A < 0.1A => No in-plane jumps'"
|
||||
print """
|
||||
|
||||
Full minimization after removal of PAFI fixes:
|
||||
energy = ${pe} eV
|
||||
|
||||
----- END PAFI ----
|
||||
"""
|
||||
|
|
@ -0,0 +1,401 @@
|
|||
LAMMPS (24 Aug 2020)
|
||||
units metal
|
||||
atom_style atomic
|
||||
atom_modify map array sort 0 0.0
|
||||
neigh_modify every 2 delay 10 check yes page 1000000 one 100000
|
||||
|
||||
## read in path data using fix property/atom- here 4th image of a NEB
|
||||
fix pa all property/atom d_nx d_ny d_nz d_dnx d_dny d_dnz d_ddnx d_ddny d_ddnz
|
||||
read_data pafipath.4.data fix pa NULL PafiPath
|
||||
Reading data file ...
|
||||
orthogonal box = (0.0000000 0.0000000 0.0000000) to (19.987184 19.987184 19.987184)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
687 atoms
|
||||
read_data CPU = 0.007 seconds
|
||||
|
||||
## EAM potential
|
||||
pair_style eam/fs
|
||||
pair_coeff * * ../../../../potentials/Fe_mm.eam.fs Fe
|
||||
Reading eam/fs potential file ../../../../potentials/Fe_mm.eam.fs with DATE: 2007-06-11
|
||||
|
||||
mass * 55.85
|
||||
thermo 100
|
||||
min_style fire
|
||||
|
||||
|
||||
compute pe all pe
|
||||
variable pe equal pe
|
||||
run 0
|
||||
Neighbor list info ...
|
||||
update every 2 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 100000, page size: 1000000
|
||||
master list distance cutoff = 7.3
|
||||
ghost atom cutoff = 7.3
|
||||
binsize = 3.65, bins = 6 6 6
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair eam/fs, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 6.903 | 6.903 | 6.903 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -2828.2492 0 -2828.2492 3015.2014
|
||||
Loop time of 6.88e-07 on 1 procs for 0 steps with 687 atoms
|
||||
|
||||
290.7% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 6.88e-07 | | |100.00
|
||||
|
||||
Nlocal: 687.000 ave 687 max 687 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 3239.00 ave 3239 max 3239 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 46843.0 ave 46843 max 46843 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 46843
|
||||
Ave neighs/atom = 68.184862
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
print "energy=${pe}"
|
||||
energy=-2828.24917967199
|
||||
|
||||
## compute property/atom to access relevant fields
|
||||
compute pa all property/atom d_nx d_ny d_nz d_dnx d_dny d_dnz d_ddnx d_ddny d_ddnz
|
||||
run 0
|
||||
Per MPI rank memory allocation (min/avg/max) = 6.903 | 6.903 | 6.903 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -2828.2492 0 -2828.2492 3015.2014
|
||||
Loop time of 9.41e-07 on 1 procs for 0 steps with 687 atoms
|
||||
|
||||
106.3% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 9.41e-07 | | |100.00
|
||||
|
||||
Nlocal: 687.000 ave 687 max 687 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 3239.00 ave 3239 max 3239 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 46843.0 ave 46843 max 46843 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 46843
|
||||
Ave neighs/atom = 68.184862
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
|
||||
## fix name group-id pafi compute-id temperature tdamp seed overdamped 0/1 com 0/1
|
||||
fix hp all pafi pa 500.0 0.01 434 overdamped no com yes
|
||||
fix pafi compute name,style: pa,property/atom
|
||||
|
||||
run 0
|
||||
Per MPI rank memory allocation (min/avg/max) = 8.403 | 8.403 | 8.403 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -2828.2492 0 -2828.2492 3015.2014
|
||||
Loop time of 5.53002e-07 on 1 procs for 0 steps with 687 atoms
|
||||
|
||||
361.7% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 5.53e-07 | | |100.00
|
||||
|
||||
Nlocal: 687.000 ave 687 max 687 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 3239.00 ave 3239 max 3239 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 46843.0 ave 46843 max 46843 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 46843
|
||||
Ave neighs/atom = 68.184862
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
minimize 0 0 1000 1000 # best if using NEB path
|
||||
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:186)
|
||||
Parameters for fire:
|
||||
dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin integrator halfstepback
|
||||
0.1 20 1.1 0.5 0.25 0.99 10 0.02 eulerimplicit yes
|
||||
Per MPI rank memory allocation (min/avg/max) = 8.403 | 8.403 | 8.403 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -2828.2492 0 -2828.2492 3015.2014
|
||||
100 9.5973821e-06 -2828.3565 0 -2828.3565 3092.936
|
||||
200 1.5010519e-13 -2828.3565 0 -2828.3565 3092.8883
|
||||
300 5.4507205e-14 -2828.3565 0 -2828.3565 3092.8883
|
||||
400 5.0283348e-14 -2828.3565 0 -2828.3565 3092.8883
|
||||
500 4.9056964e-14 -2828.3565 0 -2828.3565 3092.8883
|
||||
600 4.8645525e-14 -2828.3565 0 -2828.3565 3092.8883
|
||||
700 4.8498912e-14 -2828.3565 0 -2828.3565 3092.8882
|
||||
800 4.8444552e-14 -2828.3565 0 -2828.3565 3092.8882
|
||||
900 4.8424643e-14 -2828.3565 0 -2828.3565 3092.8882
|
||||
1000 4.84169e-14 -2828.3565 0 -2828.3565 3092.8882
|
||||
Loop time of 1.23105 on 1 procs for 1000 steps with 687 atoms
|
||||
|
||||
100.0% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
Minimization stats:
|
||||
Stopping criterion = max iterations
|
||||
Energy initial, next-to-last, final =
|
||||
-2828.24917967199 -2828.35645843665 -2828.35645844395
|
||||
Force two-norm initial, final = 1.2831004431533395 7.468917431980592e-07
|
||||
Force max component initial, final = 0.3331705829071986 2.93217567864934e-07
|
||||
Final line search alpha, max atom move = 0.0 0.0
|
||||
Iterations, force evaluations = 1000 1007
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 1.1849 | 1.1849 | 1.1849 | 0.0 | 96.25
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.010743 | 0.010743 | 0.010743 | 0.0 | 0.87
|
||||
Output | 0.00013686 | 0.00013686 | 0.00013686 | 0.0 | 0.01
|
||||
Modify | 0.02476 | 0.02476 | 0.02476 | 0.0 | 2.01
|
||||
Other | | 0.01055 | | | 0.86
|
||||
|
||||
Nlocal: 687.000 ave 687 max 687 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 3239.00 ave 3239 max 3239 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 46843.0 ave 46843 max 46843 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 46843
|
||||
Ave neighs/atom = 68.184862
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
compute dx all displace/atom
|
||||
compute dmag all reduce max c_dx[4]
|
||||
variable maxjump equal sqrt(c_dmag)
|
||||
|
||||
thermo_style custom step pe temp v_maxjump
|
||||
|
||||
|
||||
variable dn equal f_hp[4]*f_hp[4] # should be zero to float precision
|
||||
|
||||
|
||||
fix af all ave/time 1 1000 1000 f_hp[1] f_hp[2] f_hp[3] v_dn
|
||||
|
||||
variable adn equal sqrt(f_af[4]) # should be very small (approx. 1e-6 A)
|
||||
variable apf equal f_af[1]
|
||||
|
||||
|
||||
run 1000
|
||||
Per MPI rank memory allocation (min/avg/max) = 8.778 | 8.778 | 8.778 Mbytes
|
||||
Step PotEng Temp v_maxjump
|
||||
1000 -2828.3565 4.84169e-14 3.9400009e-08
|
||||
1100 -2788.9297 505.16169 0.53846275
|
||||
1200 -2784.1274 471.52773 0.5820714
|
||||
1300 -2782.7071 503.23183 0.57990998
|
||||
1400 -2780.7259 531.16141 0.61955057
|
||||
1500 -2781.512 490.9404 0.58663715
|
||||
1600 -2783.4961 513.37777 0.65062082
|
||||
1700 -2780.0312 492.76259 0.5785296
|
||||
1800 -2781.6555 512.32986 0.58277945
|
||||
1900 -2782.1648 486.16736 0.58913114
|
||||
2000 -2781.1942 493.15046 0.56962281
|
||||
Loop time of 1.3463 on 1 procs for 1000 steps with 687 atoms
|
||||
|
||||
Performance: 64.176 ns/day, 0.374 hours/ns, 742.778 timesteps/s
|
||||
100.0% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 1.2429 | 1.2429 | 1.2429 | 0.0 | 92.32
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.011341 | 0.011341 | 0.011341 | 0.0 | 0.84
|
||||
Output | 0.0001418 | 0.0001418 | 0.0001418 | 0.0 | 0.01
|
||||
Modify | 0.089655 | 0.089655 | 0.089655 | 0.0 | 6.66
|
||||
Other | | 0.002248 | | | 0.17
|
||||
|
||||
Nlocal: 687.000 ave 687 max 687 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 3239.00 ave 3239 max 3239 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 46844.0 ave 46844 max 46844 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 46844
|
||||
Ave neighs/atom = 68.186317
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
minimize 0 0 1000 1000
|
||||
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:186)
|
||||
Parameters for fire:
|
||||
dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin integrator halfstepback
|
||||
0.1 20 1.1 0.5 0.25 0.99 10 0.02 eulerimplicit yes
|
||||
Per MPI rank memory allocation (min/avg/max) = 9.278 | 9.278 | 9.278 Mbytes
|
||||
Step PotEng Temp v_maxjump
|
||||
2000 -2781.1942 0 0.56962281
|
||||
2100 -2828.3564 2.3533517e-06 0.028829312
|
||||
2200 -2828.3565 2.3166114e-09 0.0015098779
|
||||
2300 -2828.3565 5.7159208e-14 0.0017926087
|
||||
2400 -2828.3565 5.0928429e-14 0.0020264648
|
||||
2500 -2828.3565 4.9253427e-14 0.0022344754
|
||||
2600 -2828.3565 4.8706717e-14 0.0024243813
|
||||
2700 -2828.3565 4.8513392e-14 0.0026003645
|
||||
2800 -2828.3565 4.8442915e-14 0.0027651416
|
||||
2900 -2828.3565 4.8417422e-14 0.0029206277
|
||||
3000 -2828.3565 4.840734e-14 0.0030682417
|
||||
Loop time of 1.24896 on 1 procs for 1000 steps with 687 atoms
|
||||
|
||||
100.0% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
Minimization stats:
|
||||
Stopping criterion = max iterations
|
||||
Energy initial, next-to-last, final =
|
||||
-2781.19422468903 -2828.35646618005 -2828.35646618733
|
||||
Force two-norm initial, final = 29.0662462838641 7.469228601686482e-07
|
||||
Force max component initial, final = 2.2910705425983147 2.932332034966123e-07
|
||||
Final line search alpha, max atom move = 0.0 0.0
|
||||
Iterations, force evaluations = 1000 1008
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 1.2027 | 1.2027 | 1.2027 | 0.0 | 96.30
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.010706 | 0.010706 | 0.010706 | 0.0 | 0.86
|
||||
Output | 0.00016364 | 0.00016364 | 0.00016364 | 0.0 | 0.01
|
||||
Modify | 0.024873 | 0.024873 | 0.024873 | 0.0 | 1.99
|
||||
Other | | 0.01053 | | | 0.84
|
||||
|
||||
Nlocal: 687.000 ave 687 max 687 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 3209.00 ave 3209 max 3209 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 48576.0 ave 48576 max 48576 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 48576
|
||||
Ave neighs/atom = 70.707424
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
|
||||
variable s_pe equal ${pe}
|
||||
variable s_pe equal -2828.35646618733
|
||||
variable s_apf equal ${apf}
|
||||
variable s_apf equal -0.117797263023891
|
||||
variable s_adn equal ${adn}
|
||||
variable s_adn equal 1.36645874805011e-05
|
||||
variable s_maxjump equal ${maxjump}
|
||||
variable s_maxjump equal 0.00306824174470598
|
||||
|
||||
unfix hp
|
||||
unfix pa
|
||||
unfix af
|
||||
minimize 0 0 1000 1000
|
||||
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:186)
|
||||
Parameters for fire:
|
||||
dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin integrator halfstepback
|
||||
0.1 20 1.1 0.5 0.25 0.99 10 0.02 eulerimplicit yes
|
||||
Per MPI rank memory allocation (min/avg/max) = 8.778 | 8.778 | 8.778 Mbytes
|
||||
Step PotEng Temp v_maxjump
|
||||
3000 -2828.3565 0 0.0030682417
|
||||
3100 -2828.5793 0.00016081731 0.54349309
|
||||
3200 -2828.5793 1.1620454e-14 0.54386121
|
||||
3300 -2828.5793 4.4699993e-19 0.54386124
|
||||
3400 -2828.5793 1.0681196e-23 0.54386124
|
||||
3500 -2828.5793 2.046642e-25 0.54386124
|
||||
3600 -2828.5793 1.4539039e-25 0.54386124
|
||||
3700 -2828.5793 2.0152508e-25 0.54386124
|
||||
3800 -2828.5793 1.3593174e-25 0.54386124
|
||||
3900 -2828.5793 1.8155665e-25 0.54386124
|
||||
4000 -2828.5793 1.2308872e-25 0.54386124
|
||||
Loop time of 1.20543 on 1 procs for 1000 steps with 687 atoms
|
||||
|
||||
100.0% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
Minimization stats:
|
||||
Stopping criterion = max iterations
|
||||
Energy initial, next-to-last, final =
|
||||
-2828.35646618733 -2828.57925588226 -2828.57925588226
|
||||
Force two-norm initial, final = 0.5639141583763609 2.3523921777846766e-13
|
||||
Force max component initial, final = 0.22089321768110187 2.037953139577553e-14
|
||||
Final line search alpha, max atom move = 0.0 0.0
|
||||
Iterations, force evaluations = 1000 1017
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 1.1843 | 1.1843 | 1.1843 | 0.0 | 98.25
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.010469 | 0.010469 | 0.010469 | 0.0 | 0.87
|
||||
Output | 0.00012213 | 0.00012213 | 0.00012213 | 0.0 | 0.01
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 0.01051 | | | 0.87
|
||||
|
||||
Nlocal: 687.000 ave 687 max 687 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 3239.00 ave 3239 max 3239 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 46844.0 ave 46844 max 46844 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 46844
|
||||
Ave neighs/atom = 68.186317
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
|
||||
|
||||
print """
|
||||
---------- PAFI RESULTS --------
|
||||
During run:
|
||||
Average Distance From Hyperplane = ${s_adn}A (Should be very small, around 1e-5A),
|
||||
Average Force Projection Along MFEP = ${s_apf} eV/A
|
||||
"""
|
||||
|
||||
---------- PAFI RESULTS --------
|
||||
During run:
|
||||
Average Distance From Hyperplane = 1.36645874805011e-05A (Should be very small, around 1e-5A),
|
||||
Average Force Projection Along MFEP = -0.117797263023891 eV/A
|
||||
|
||||
print """
|
||||
In-plane minimization post-run:
|
||||
energy = ${s_pe} eV
|
||||
"""
|
||||
|
||||
In-plane minimization post-run:
|
||||
energy = -2828.35646618733 eV
|
||||
|
||||
if "${s_maxjump} > 0.1" then "print 'Max Atomic Displacement = ${s_maxjump}A > 0.1A'" "print ' => Possible shallow in-plane metastability: Reduce time in basin and/or decrease tdamp'" else "print 'Max Atomic Displacement = ${s_maxjump}A < 0.1A => No in-plane jumps'"
|
||||
print 'Max Atomic Displacement = ${s_maxjump}A < 0.1A => No in-plane jumps'
|
||||
Max Atomic Displacement = 0.00306824174470598A < 0.1A => No in-plane jumps
|
||||
print """
|
||||
Full minimization after removal of PAFI fixes:
|
||||
energy = ${pe} eV
|
||||
----- END PAFI ----
|
||||
"""
|
||||
|
||||
Full minimization after removal of PAFI fixes:
|
||||
energy = -2828.57925588226 eV
|
||||
----- END PAFI ----
|
||||
|
||||
|
||||
Please see the log.cite file for references relevant to this simulation
|
||||
|
||||
Total wall time: 0:00:05
|
||||
|
|
@ -0,0 +1,401 @@
|
|||
LAMMPS (24 Aug 2020)
|
||||
units metal
|
||||
atom_style atomic
|
||||
atom_modify map array sort 0 0.0
|
||||
neigh_modify every 2 delay 10 check yes page 1000000 one 100000
|
||||
|
||||
## read in path data using fix property/atom- here 4th image of a NEB
|
||||
fix pa all property/atom d_nx d_ny d_nz d_dnx d_dny d_dnz d_ddnx d_ddny d_ddnz
|
||||
read_data pafipath.4.data fix pa NULL PafiPath
|
||||
Reading data file ...
|
||||
orthogonal box = (0.0000000 0.0000000 0.0000000) to (19.987184 19.987184 19.987184)
|
||||
2 by 1 by 2 MPI processor grid
|
||||
reading atoms ...
|
||||
687 atoms
|
||||
read_data CPU = 0.007 seconds
|
||||
|
||||
## EAM potential
|
||||
pair_style eam/fs
|
||||
pair_coeff * * ../../../../potentials/Fe_mm.eam.fs Fe
|
||||
Reading eam/fs potential file ../../../../potentials/Fe_mm.eam.fs with DATE: 2007-06-11
|
||||
|
||||
mass * 55.85
|
||||
thermo 100
|
||||
min_style fire
|
||||
|
||||
|
||||
compute pe all pe
|
||||
variable pe equal pe
|
||||
run 0
|
||||
Neighbor list info ...
|
||||
update every 2 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 100000, page size: 1000000
|
||||
master list distance cutoff = 7.3
|
||||
ghost atom cutoff = 7.3
|
||||
binsize = 3.65, bins = 6 6 6
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair eam/fs, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 6.844 | 6.844 | 6.844 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -2828.2492 0 -2828.2492 3015.2014
|
||||
Loop time of 2.82e-06 on 4 procs for 0 steps with 687 atoms
|
||||
|
||||
97.5% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 2.82e-06 | | |100.00
|
||||
|
||||
Nlocal: 171.750 ave 187 max 163 min
|
||||
Histogram: 1 1 1 0 0 0 0 0 0 1
|
||||
Nghost: 1854.25 ave 1863 max 1839 min
|
||||
Histogram: 1 0 0 0 0 0 0 1 1 1
|
||||
Neighs: 11710.8 ave 12319 max 10647 min
|
||||
Histogram: 1 0 0 0 0 0 0 1 1 1
|
||||
|
||||
Total # of neighbors = 46843
|
||||
Ave neighs/atom = 68.184862
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
print "energy=${pe}"
|
||||
energy=-2828.24917967201
|
||||
|
||||
## compute property/atom to access relevant fields
|
||||
compute pa all property/atom d_nx d_ny d_nz d_dnx d_dny d_dnz d_ddnx d_ddny d_ddnz
|
||||
run 0
|
||||
Per MPI rank memory allocation (min/avg/max) = 6.844 | 6.844 | 6.844 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -2828.2492 0 -2828.2492 3015.2014
|
||||
Loop time of 2.667e-06 on 4 procs for 0 steps with 687 atoms
|
||||
|
||||
37.5% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 2.667e-06 | | |100.00
|
||||
|
||||
Nlocal: 171.750 ave 187 max 163 min
|
||||
Histogram: 1 1 1 0 0 0 0 0 0 1
|
||||
Nghost: 1854.25 ave 1863 max 1839 min
|
||||
Histogram: 1 0 0 0 0 0 0 1 1 1
|
||||
Neighs: 11710.8 ave 12319 max 10647 min
|
||||
Histogram: 1 0 0 0 0 0 0 1 1 1
|
||||
|
||||
Total # of neighbors = 46843
|
||||
Ave neighs/atom = 68.184862
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
|
||||
## fix name group-id pafi compute-id temperature tdamp seed overdamped 0/1 com 0/1
|
||||
fix hp all pafi pa 500.0 0.01 434 overdamped no com yes
|
||||
fix pafi compute name,style: pa,property/atom
|
||||
|
||||
run 0
|
||||
Per MPI rank memory allocation (min/avg/max) = 8.344 | 8.344 | 8.344 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -2828.2492 0 -2828.2492 3015.2014
|
||||
Loop time of 2.52725e-06 on 4 procs for 0 steps with 687 atoms
|
||||
|
||||
108.8% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 2.527e-06 | | |100.00
|
||||
|
||||
Nlocal: 171.750 ave 187 max 163 min
|
||||
Histogram: 1 1 1 0 0 0 0 0 0 1
|
||||
Nghost: 1854.25 ave 1863 max 1839 min
|
||||
Histogram: 1 0 0 0 0 0 0 1 1 1
|
||||
Neighs: 11710.8 ave 12319 max 10647 min
|
||||
Histogram: 1 0 0 0 0 0 0 1 1 1
|
||||
|
||||
Total # of neighbors = 46843
|
||||
Ave neighs/atom = 68.184862
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
minimize 0 0 1000 1000 # best if using NEB path
|
||||
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:186)
|
||||
Parameters for fire:
|
||||
dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin integrator halfstepback
|
||||
0.1 20 1.1 0.5 0.25 0.99 10 0.02 eulerimplicit yes
|
||||
Per MPI rank memory allocation (min/avg/max) = 8.344 | 8.344 | 8.344 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -2828.2492 0 -2828.2492 3015.2014
|
||||
100 9.5973821e-06 -2828.3565 0 -2828.3565 3092.936
|
||||
200 1.501052e-13 -2828.3565 0 -2828.3565 3092.8883
|
||||
300 5.4507205e-14 -2828.3565 0 -2828.3565 3092.8883
|
||||
400 5.0283344e-14 -2828.3565 0 -2828.3565 3092.8883
|
||||
500 4.9056968e-14 -2828.3565 0 -2828.3565 3092.8883
|
||||
600 4.8645528e-14 -2828.3565 0 -2828.3565 3092.8883
|
||||
700 4.8498892e-14 -2828.3565 0 -2828.3565 3092.8882
|
||||
800 4.8444568e-14 -2828.3565 0 -2828.3565 3092.8882
|
||||
900 4.8424653e-14 -2828.3565 0 -2828.3565 3092.8882
|
||||
1000 4.8416878e-14 -2828.3565 0 -2828.3565 3092.8882
|
||||
Loop time of 0.43333 on 4 procs for 1000 steps with 687 atoms
|
||||
|
||||
95.4% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
Minimization stats:
|
||||
Stopping criterion = max iterations
|
||||
Energy initial, next-to-last, final =
|
||||
-2828.24917967201 -2828.35645843663 -2828.3564584439
|
||||
Force two-norm initial, final = 1.2831004431533406 7.468917439994878e-07
|
||||
Force max component initial, final = 0.33317058290719853 2.9321733663323357e-07
|
||||
Final line search alpha, max atom move = 0.0 0.0
|
||||
Iterations, force evaluations = 1000 1007
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.3303 | 0.34246 | 0.35813 | 1.7 | 79.03
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.046971 | 0.062117 | 0.073961 | 3.9 | 14.33
|
||||
Output | 0.00013746 | 0.00014874 | 0.00016997 | 0.0 | 0.03
|
||||
Modify | 0.0099064 | 0.01077 | 0.011612 | 0.6 | 2.49
|
||||
Other | | 0.01784 | | | 4.12
|
||||
|
||||
Nlocal: 171.750 ave 187 max 163 min
|
||||
Histogram: 1 1 1 0 0 0 0 0 0 1
|
||||
Nghost: 1854.25 ave 1863 max 1839 min
|
||||
Histogram: 1 0 0 0 0 0 0 1 1 1
|
||||
Neighs: 11710.8 ave 12319 max 10647 min
|
||||
Histogram: 1 0 0 0 0 0 0 1 1 1
|
||||
|
||||
Total # of neighbors = 46843
|
||||
Ave neighs/atom = 68.184862
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
compute dx all displace/atom
|
||||
compute dmag all reduce max c_dx[4]
|
||||
variable maxjump equal sqrt(c_dmag)
|
||||
|
||||
thermo_style custom step pe temp v_maxjump
|
||||
|
||||
|
||||
variable dn equal f_hp[4]*f_hp[4] # should be zero to float precision
|
||||
|
||||
|
||||
fix af all ave/time 1 1000 1000 f_hp[1] f_hp[2] f_hp[3] v_dn
|
||||
|
||||
variable adn equal sqrt(f_af[4]) # should be very small (approx. 1e-6 A)
|
||||
variable apf equal f_af[1]
|
||||
|
||||
|
||||
run 1000
|
||||
Per MPI rank memory allocation (min/avg/max) = 8.719 | 8.719 | 8.719 Mbytes
|
||||
Step PotEng Temp v_maxjump
|
||||
1000 -2828.3565 4.8416878e-14 3.9625012e-08
|
||||
1100 -2787.3877 503.1749 0.53858794
|
||||
1200 -2783.8463 508.91879 0.57609619
|
||||
1300 -2782.5288 533.03961 0.64090828
|
||||
1400 -2782.4551 508.05688 0.61706248
|
||||
1500 -2784.0499 538.34202 0.56379358
|
||||
1600 -2786.4068 487.9633 0.57600511
|
||||
1700 -2780.7727 489.2309 0.61048418
|
||||
1800 -2780.4909 519.01911 0.60204913
|
||||
1900 -2782.0825 512.3441 0.6071772
|
||||
2000 -2779.5449 526.21838 0.58247597
|
||||
Loop time of 0.446418 on 4 procs for 1000 steps with 687 atoms
|
||||
|
||||
Performance: 193.541 ns/day, 0.124 hours/ns, 2240.053 timesteps/s
|
||||
95.8% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.34383 | 0.35018 | 0.35655 | 0.9 | 78.44
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.049244 | 0.055398 | 0.061404 | 2.2 | 12.41
|
||||
Output | 0.00015376 | 0.00016326 | 0.00018954 | 0.0 | 0.04
|
||||
Modify | 0.034111 | 0.034443 | 0.03473 | 0.1 | 7.72
|
||||
Other | | 0.006229 | | | 1.40
|
||||
|
||||
Nlocal: 171.750 ave 183 max 165 min
|
||||
Histogram: 1 1 0 1 0 0 0 0 0 1
|
||||
Nghost: 1854.25 ave 1861 max 1843 min
|
||||
Histogram: 1 0 0 0 0 0 1 0 1 1
|
||||
Neighs: 11711.0 ave 12115 max 10983 min
|
||||
Histogram: 1 0 0 0 0 0 1 0 1 1
|
||||
|
||||
Total # of neighbors = 46844
|
||||
Ave neighs/atom = 68.186317
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
minimize 0 0 1000 1000
|
||||
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:186)
|
||||
Parameters for fire:
|
||||
dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin integrator halfstepback
|
||||
0.1 20 1.1 0.5 0.25 0.99 10 0.02 eulerimplicit yes
|
||||
Per MPI rank memory allocation (min/avg/max) = 9.219 | 9.219 | 9.219 Mbytes
|
||||
Step PotEng Temp v_maxjump
|
||||
2000 -2779.5449 0 0.58247597
|
||||
2100 -2828.3564 5.5901288e-06 0.022258122
|
||||
2200 -2828.3565 2.5724945e-11 0.0014462408
|
||||
2300 -2828.3565 5.8980018e-14 0.0018217607
|
||||
2400 -2828.3565 5.1597335e-14 0.0020531463
|
||||
2500 -2828.3565 4.9451623e-14 0.0022588465
|
||||
2600 -2828.3565 4.8775183e-14 0.0024468988
|
||||
2700 -2828.3565 4.8538707e-14 0.0026213819
|
||||
2800 -2828.3565 4.8450209e-14 0.0027849194
|
||||
2900 -2828.3565 4.8421002e-14 0.0029393609
|
||||
3000 -2828.3565 4.8409037e-14 0.0030860795
|
||||
Loop time of 0.432206 on 4 procs for 1000 steps with 687 atoms
|
||||
|
||||
91.9% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
Minimization stats:
|
||||
Stopping criterion = max iterations
|
||||
Energy initial, next-to-last, final =
|
||||
-2779.54493572831 -2828.35646627033 -2828.3564662776
|
||||
Force two-norm initial, final = 29.30132707276523 7.469232246845132e-07
|
||||
Force max component initial, final = 2.0790292865239595 2.9323707292916446e-07
|
||||
Final line search alpha, max atom move = 0.0 0.0
|
||||
Iterations, force evaluations = 1000 1008
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.33828 | 0.34752 | 0.35566 | 1.1 | 80.41
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.04768 | 0.056242 | 0.065447 | 2.8 | 13.01
|
||||
Output | 0.00016892 | 0.00018009 | 0.00020427 | 0.0 | 0.04
|
||||
Modify | 0.010016 | 0.010698 | 0.011316 | 0.4 | 2.48
|
||||
Other | | 0.01756 | | | 4.06
|
||||
|
||||
Nlocal: 171.750 ave 184 max 165 min
|
||||
Histogram: 1 0 2 0 0 0 0 0 0 1
|
||||
Nghost: 1828.25 ave 1835 max 1815 min
|
||||
Histogram: 1 0 0 0 0 0 0 1 0 2
|
||||
Neighs: 12157.5 ave 12578 max 11155 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 3
|
||||
|
||||
Total # of neighbors = 48630
|
||||
Ave neighs/atom = 70.786026
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
|
||||
variable s_pe equal ${pe}
|
||||
variable s_pe equal -2828.3564662776
|
||||
variable s_apf equal ${apf}
|
||||
variable s_apf equal -0.44222069537824
|
||||
variable s_adn equal ${adn}
|
||||
variable s_adn equal 1.36645875188617e-05
|
||||
variable s_maxjump equal ${maxjump}
|
||||
variable s_maxjump equal 0.00308607949485962
|
||||
|
||||
unfix hp
|
||||
unfix pa
|
||||
unfix af
|
||||
minimize 0 0 1000 1000
|
||||
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:186)
|
||||
Parameters for fire:
|
||||
dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin integrator halfstepback
|
||||
0.1 20 1.1 0.5 0.25 0.99 10 0.02 eulerimplicit yes
|
||||
Per MPI rank memory allocation (min/avg/max) = 8.719 | 8.719 | 8.719 Mbytes
|
||||
Step PotEng Temp v_maxjump
|
||||
3000 -2828.3565 0 0.0030860795
|
||||
3100 -2828.5793 0.00016081696 0.5434931
|
||||
3200 -2828.5793 1.1620338e-14 0.54386121
|
||||
3300 -2828.5793 4.4708247e-19 0.54386124
|
||||
3400 -2828.5793 1.0642423e-23 0.54386124
|
||||
3500 -2828.5793 1.9858026e-25 0.54386124
|
||||
3600 -2828.5793 1.395464e-25 0.54386124
|
||||
3700 -2828.5793 2.0687121e-25 0.54386124
|
||||
3800 -2828.5793 1.3586688e-25 0.54386124
|
||||
3900 -2828.5793 2.0485828e-25 0.54386124
|
||||
4000 -2828.5793 1.3803592e-25 0.54386124
|
||||
Loop time of 0.427843 on 4 procs for 1000 steps with 687 atoms
|
||||
|
||||
96.5% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
Minimization stats:
|
||||
Stopping criterion = max iterations
|
||||
Energy initial, next-to-last, final =
|
||||
-2828.3564662776 -2828.57925588218 -2828.57925588218
|
||||
Force two-norm initial, final = 0.5639144322930092 2.34573110665124e-13
|
||||
Force max component initial, final = 0.22089332498188619 2.0969337377607644e-14
|
||||
Final line search alpha, max atom move = 0.0 0.0
|
||||
Iterations, force evaluations = 1000 1017
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.34748 | 0.35476 | 0.3611 | 0.9 | 82.92
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.049149 | 0.054906 | 0.062132 | 2.2 | 12.83
|
||||
Output | 0.00015435 | 0.00016534 | 0.00018687 | 0.0 | 0.04
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 0.01801 | | | 4.21
|
||||
|
||||
Nlocal: 171.750 ave 183 max 165 min
|
||||
Histogram: 1 1 0 1 0 0 0 0 0 1
|
||||
Nghost: 1854.25 ave 1861 max 1843 min
|
||||
Histogram: 1 0 0 0 0 0 1 0 1 1
|
||||
Neighs: 11711.0 ave 12115 max 10983 min
|
||||
Histogram: 1 0 0 0 0 0 1 0 1 1
|
||||
|
||||
Total # of neighbors = 46844
|
||||
Ave neighs/atom = 68.186317
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
|
||||
|
||||
print """
|
||||
---------- PAFI RESULTS --------
|
||||
During run:
|
||||
Average Distance From Hyperplane = ${s_adn}A (Should be very small, around 1e-5A),
|
||||
Average Force Projection Along MFEP = ${s_apf} eV/A
|
||||
"""
|
||||
|
||||
---------- PAFI RESULTS --------
|
||||
During run:
|
||||
Average Distance From Hyperplane = 1.36645875188617e-05A (Should be very small, around 1e-5A),
|
||||
Average Force Projection Along MFEP = -0.44222069537824 eV/A
|
||||
|
||||
print """
|
||||
In-plane minimization post-run:
|
||||
energy = ${s_pe} eV
|
||||
"""
|
||||
|
||||
In-plane minimization post-run:
|
||||
energy = -2828.3564662776 eV
|
||||
|
||||
if "${s_maxjump} > 0.1" then "print 'Max Atomic Displacement = ${s_maxjump}A > 0.1A'" "print ' => Possible shallow in-plane metastability: Reduce time in basin and/or decrease tdamp'" else "print 'Max Atomic Displacement = ${s_maxjump}A < 0.1A => No in-plane jumps'"
|
||||
print 'Max Atomic Displacement = ${s_maxjump}A < 0.1A => No in-plane jumps'
|
||||
Max Atomic Displacement = 0.00308607949485962A < 0.1A => No in-plane jumps
|
||||
print """
|
||||
Full minimization after removal of PAFI fixes:
|
||||
energy = ${pe} eV
|
||||
----- END PAFI ----
|
||||
"""
|
||||
|
||||
Full minimization after removal of PAFI fixes:
|
||||
energy = -2828.57925588218 eV
|
||||
----- END PAFI ----
|
||||
|
||||
|
||||
Please see the log.cite file for references relevant to this simulation
|
||||
|
||||
Total wall time: 0:00:01
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -270,6 +270,27 @@ class lammps(object):
|
|||
[c_void_p,c_char_p,c_int,c_int,c_int,POINTER(c_int),c_void_p]
|
||||
self.lib.lammps_scatter_atoms_subset.restype = None
|
||||
|
||||
self.lib.lammps_gather.argtypes = \
|
||||
[c_void_p,c_char_p,c_int,c_int,c_void_p]
|
||||
self.lib.lammps_gather.restype = None
|
||||
|
||||
self.lib.lammps_gather_concat.argtypes = \
|
||||
[c_void_p,c_char_p,c_int,c_int,c_void_p]
|
||||
self.lib.lammps_gather_concat.restype = None
|
||||
|
||||
self.lib.lammps_gather_subset.argtypes = \
|
||||
[c_void_p,c_char_p,c_int,c_int,c_int,POINTER(c_int),c_void_p]
|
||||
self.lib.lammps_gather_subset.restype = None
|
||||
|
||||
self.lib.lammps_scatter.argtypes = \
|
||||
[c_void_p,c_char_p,c_int,c_int,c_void_p]
|
||||
self.lib.lammps_scatter.restype = None
|
||||
|
||||
self.lib.lammps_scatter_subset.argtypes = \
|
||||
[c_void_p,c_char_p,c_int,c_int,c_int,POINTER(c_int),c_void_p]
|
||||
self.lib.lammps_scatter_subset.restype = None
|
||||
|
||||
|
||||
self.lib.lammps_find_pair_neighlist.argtypes = [c_void_p, c_char_p, c_int, c_int, c_int]
|
||||
self.lib.lammps_find_pair_neighlist.restype = c_int
|
||||
|
||||
|
|
@ -1151,8 +1172,67 @@ class lammps(object):
|
|||
if name: name = name.encode()
|
||||
self.lib.lammps_scatter_atoms_subset(self.lmp,name,type,count,ndata,ids,data)
|
||||
|
||||
# return vector of atom/compute/fix properties gathered across procs
|
||||
# 3 variants to match src/library.cpp
|
||||
# name = atom property recognized by LAMMPS in atom->extract()
|
||||
# type = 0 for integer values, 1 for double values
|
||||
# count = number of per-atom valus, 1 for type or charge, 3 for x or f
|
||||
# returned data is a 1d vector - doc how it is ordered?
|
||||
# NOTE: need to insure are converting to/from correct Python type
|
||||
# e.g. for Python list or NumPy or ctypes
|
||||
def gather(self,name,type,count):
|
||||
if name: name = name.encode()
|
||||
natoms = self.lib.lammps_get_natoms(self.lmp)
|
||||
if type == 0:
|
||||
data = ((count*natoms)*c_int)()
|
||||
self.lib.lammps_gather(self.lmp,name,type,count,data)
|
||||
elif type == 1:
|
||||
data = ((count*natoms)*c_double)()
|
||||
self.lib.lammps_gather(self.lmp,name,type,count,data)
|
||||
else: return None
|
||||
return data
|
||||
|
||||
# -------------------------------------------------------------------------
|
||||
def gather_concat(self,name,type,count):
|
||||
if name: name = name.encode()
|
||||
natoms = self.lib.lammps_get_natoms(self.lmp)
|
||||
if type == 0:
|
||||
data = ((count*natoms)*c_int)()
|
||||
self.lib.lammps_gather_concat(self.lmp,name,type,count,data)
|
||||
elif type == 1:
|
||||
data = ((count*natoms)*c_double)()
|
||||
self.lib.lammps_gather_concat(self.lmp,name,type,count,data)
|
||||
else: return None
|
||||
return data
|
||||
|
||||
def gather_subset(self,name,type,count,ndata,ids):
|
||||
if name: name = name.encode()
|
||||
if type == 0:
|
||||
data = ((count*ndata)*c_int)()
|
||||
self.lib.lammps_gather_subset(self.lmp,name,type,count,ndata,ids,data)
|
||||
elif type == 1:
|
||||
data = ((count*ndata)*c_double)()
|
||||
self.lib.lammps_gather_subset(self.lmp,name,type,count,ndata,ids,data)
|
||||
else: return None
|
||||
return data
|
||||
|
||||
# scatter vector of atom/compute/fix properties across procs
|
||||
# 2 variants to match src/library.cpp
|
||||
# name = atom property recognized by LAMMPS in atom->extract()
|
||||
# type = 0 for integer values, 1 for double values
|
||||
# count = number of per-atom valus, 1 for type or charge, 3 for x or f
|
||||
# assume data is of correct type and length, as created by gather_atoms()
|
||||
# NOTE: need to insure are converting to/from correct Python type
|
||||
# e.g. for Python list or NumPy or ctypes
|
||||
|
||||
def scatter(self,name,type,count,data):
|
||||
if name: name = name.encode()
|
||||
self.lib.lammps_scatter(self.lmp,name,type,count,data)
|
||||
|
||||
def scatter_subset(self,name,type,count,ndata,ids,data):
|
||||
if name: name = name.encode()
|
||||
self.lib.lammps_scatter_subset(self.lmp,name,type,count,ndata,ids,data)
|
||||
|
||||
# -------------------------------------------------------------------------
|
||||
|
||||
def encode_image_flags(self,ix,iy,iz):
|
||||
""" convert 3 integers with image flags for x-, y-, and z-direction
|
||||
|
|
@ -1192,6 +1272,7 @@ class lammps(object):
|
|||
return [int(i) for i in flags]
|
||||
|
||||
# -------------------------------------------------------------------------
|
||||
|
||||
# create N atoms on all procs
|
||||
# N = global number of atoms
|
||||
# id = ID of each atom (optional, can be None)
|
||||
|
|
|
|||
|
|
@ -70,6 +70,9 @@
|
|||
/fix_*cauchy.cpp
|
||||
/fix_*cauchy.h
|
||||
|
||||
/fix_pafi*.cpp
|
||||
/fix_pafi*.h
|
||||
|
||||
/compute_test_nbl.cpp
|
||||
/compute_test_nbl.h
|
||||
/pair_multi_lucy.cpp
|
||||
|
|
|
|||
|
|
@ -63,6 +63,7 @@ fix momentum/chunk, Jiang Xiao (Hong Kong Polytechnic University), polyu-xiao.ji
|
|||
fix npt/cauchy, R. E. Miller (Carleton University), F. Pavia and S. Pattamatta, 12 Jan 2020
|
||||
fix nvk, Efrem Braun (UC Berkeley), efrem.braun at gmail.com, https://github.com/lammps/lammps/pull/310
|
||||
fix orient/eco Adrian A. Schratt and Volker Mohles (Ruhr-Uni Bochum), volker.mohles at rub.de, 6 Jun 2020
|
||||
fix pafi, Thomas Swinburne (CNRS), swinburne at cinam.univ-mrs.fr, 1st Sep 2020
|
||||
fix pimd, Yuxing Peng (U Chicago), yuxing at uchicago.edu, 24 Nov 2014
|
||||
fix propel/self, Stefan Paquay (Brandeis U), stefanpaquay at gmail.com, 20 Jan 2020
|
||||
fix rhok, Ulf Pedersen (Roskilde U), ulf at urp.dk, 25 Sep 2017
|
||||
|
|
|
|||
|
|
@ -0,0 +1,751 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 authors: Thomas Swinburne (CNRS & CINaM, Marseille, France)
|
||||
|
||||
Please cite the related publication:
|
||||
T.D. Swinburne and M.-C. Marinica, Unsupervised calculation of free energy barriers in large crystalline systems, Physical Review Letters 2018
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
|
||||
#include "fix_pafi.h"
|
||||
#include <mpi.h>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <cstdlib>
|
||||
#include "math_extra.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "update.h"
|
||||
#include "modify.h"
|
||||
#include "domain.h"
|
||||
#include "respa.h"
|
||||
#include "comm.h"
|
||||
#include "compute.h"
|
||||
#include "random_mars.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "utils.h"
|
||||
#include "citeme.h"
|
||||
#include "fmt/format.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
static const char cite_fix_pafi_package[] =
|
||||
"citation for fix pafi:\n\n"
|
||||
"@article{SwinburneMarinica2018,\n"
|
||||
"author={T. D. Swinburne and M. C. Marinica},\n"
|
||||
"title={Unsupervised calculation of free energy barriers in large "
|
||||
"crystalline systems},\n"
|
||||
"journal={Physical Review Letters},\n"
|
||||
"volume={120},\n"
|
||||
"number={13},\n"
|
||||
"pages={135503},\n"
|
||||
"year={2018},\n"
|
||||
"publisher={APS}\n}\n"
|
||||
"Recommended to be coupled with PAFI++ code:\n"
|
||||
"https://github.com/tomswinburne/pafi\n";
|
||||
|
||||
using namespace FixConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixPAFI::FixPAFI(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg), random(NULL), computename(NULL),
|
||||
h(NULL), step_respa(NULL)
|
||||
{
|
||||
if (lmp->citeme) lmp->citeme->add(cite_fix_pafi_package);
|
||||
|
||||
if (narg < 11) error->all(FLERR,"Illegal fix pafi command");
|
||||
|
||||
dynamic_group_allow = 0;
|
||||
vector_flag = 1;
|
||||
size_vector = 5;
|
||||
global_freq = 1;
|
||||
extvector = 0;
|
||||
od_flag = 0;
|
||||
com_flag = 1;
|
||||
time_integrate = 1;
|
||||
|
||||
int n = strlen(arg[3])+1;
|
||||
computename = new char[n];
|
||||
strcpy(computename,&arg[3][0]);
|
||||
|
||||
icompute = modify->find_compute(computename);
|
||||
if (icompute == -1)
|
||||
error->all(FLERR,"Compute ID for fix pafi does not exist");
|
||||
PathCompute = modify->compute[icompute];
|
||||
if (PathCompute->peratom_flag==0)
|
||||
error->all(FLERR,"Compute for fix pafi does not calculate a local array");
|
||||
if (PathCompute->size_peratom_cols < 9)
|
||||
error->all(FLERR,"Compute for fix pafi must have 9 fields per atom");
|
||||
|
||||
if (comm->me==0)
|
||||
utils::logmesg(lmp,fmt::format("fix pafi compute name,style: {},{}\n",computename,PathCompute->style));
|
||||
|
||||
respa_level_support = 1;
|
||||
ilevel_respa = nlevels_respa = 0;
|
||||
|
||||
temperature = utils::numeric(FLERR,arg[4],false,lmp);
|
||||
t_period = utils::numeric(FLERR,arg[5],false,lmp);
|
||||
seed = utils::inumeric(FLERR,arg[6],false,lmp);
|
||||
// TODO UNITS
|
||||
gamma = 1. / t_period / force->ftm2v;
|
||||
sqrtD = sqrt(1.) * sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e*temperature) / force->ftm2v;
|
||||
|
||||
// optional args
|
||||
int iarg = 7;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"overdamped") == 0) {
|
||||
if (strcmp(arg[iarg+1],"no") == 0) od_flag = 0;
|
||||
else if (strcmp(arg[iarg+1],"0") == 0) od_flag = 0;
|
||||
else if (strcmp(arg[iarg+1],"yes") == 0) od_flag = 1;
|
||||
else if (strcmp(arg[iarg+1],"1") == 0) od_flag = 1;
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"com") == 0) {
|
||||
if (strcmp(arg[iarg+1],"no") == 0) com_flag = 0;
|
||||
else if (strcmp(arg[iarg+1],"0") == 0) com_flag = 0;
|
||||
else if (strcmp(arg[iarg+1],"yes") == 0) com_flag = 1;
|
||||
else if (strcmp(arg[iarg+1],"1") == 0) com_flag = 1;
|
||||
iarg += 2;
|
||||
} else error->all(FLERR,"Illegal fix pafi command");
|
||||
}
|
||||
force_flag = 0;
|
||||
|
||||
for(int i = 0; i < 10; i++) {
|
||||
c_v[i] = 0.;
|
||||
c_v_all[i] = 0.;
|
||||
}
|
||||
for(int i=0; i<6; i++) {
|
||||
proj[i] = 0.0;
|
||||
proj_all[i] = 0.0;
|
||||
}
|
||||
for(int i=0; i<5; i++) {
|
||||
results[i] = 0.0;
|
||||
results_all[i] = 0.0;
|
||||
}
|
||||
maxatom = 1;
|
||||
memory->create(h,maxatom,3,"FixPAFI:h");
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
random = new RanMars(lmp,seed + comm->me);
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixPAFI::~FixPAFI()
|
||||
{
|
||||
if (copymode) return;
|
||||
delete random;
|
||||
delete [] computename;
|
||||
memory->destroy(h);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixPAFI::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= POST_FORCE;
|
||||
mask |= POST_FORCE_RESPA;
|
||||
mask |= MIN_POST_FORCE;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
// nve
|
||||
mask |= FINAL_INTEGRATE;
|
||||
mask |= INITIAL_INTEGRATE_RESPA;
|
||||
mask |= FINAL_INTEGRATE_RESPA;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixPAFI::init()
|
||||
{
|
||||
dtv = update->dt;
|
||||
dtf = 0.5 * update->dt * force->ftm2v;
|
||||
|
||||
icompute = modify->find_compute(computename);
|
||||
if (icompute == -1)
|
||||
error->all(FLERR,"Compute ID for fix pafi does not exist");
|
||||
PathCompute = modify->compute[icompute];
|
||||
if (PathCompute->peratom_flag==0)
|
||||
error->all(FLERR,"Compute for fix pafi does not calculate a local array");
|
||||
if (PathCompute->size_peratom_cols < 9)
|
||||
error->all(FLERR,"Compute for fix pafi must have 9 fields per atom");
|
||||
|
||||
|
||||
if (strstr(update->integrate_style,"respa")) {
|
||||
step_respa = ((Respa *) update->integrate)->step; // nve
|
||||
nlevels_respa = ((Respa *) update->integrate)->nlevels;
|
||||
if (respa_level >= 0) ilevel_respa = MIN(respa_level,nlevels_respa-1);
|
||||
else ilevel_respa = nlevels_respa-1;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void FixPAFI::setup(int vflag)
|
||||
{
|
||||
if (strstr(update->integrate_style,"verlet"))
|
||||
post_force(vflag);
|
||||
else
|
||||
for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
|
||||
((Respa *) update->integrate)->copy_flevel_f(ilevel);
|
||||
post_force_respa(vflag,ilevel,0);
|
||||
((Respa *) update->integrate)->copy_f_flevel(ilevel);
|
||||
}
|
||||
}
|
||||
|
||||
void FixPAFI::min_setup(int vflag)
|
||||
{
|
||||
if( utils::strmatch(update->minimize_style,"^fire")==0 &&
|
||||
utils::strmatch(update->minimize_style,"^quickmin")==0 )
|
||||
error->all(FLERR,"fix pafi requires a damped dynamics minimizer");
|
||||
min_post_force(vflag);
|
||||
}
|
||||
|
||||
|
||||
void FixPAFI::post_force(int vflag)
|
||||
{
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double *rmass = atom->rmass;
|
||||
double *mass = atom->mass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
// reallocate norm array if necessary
|
||||
if (atom->nmax > maxatom) {
|
||||
maxatom = atom->nmax;
|
||||
memory->destroy(h);
|
||||
memory->create(h,maxatom,3,"FixPAFI:h");
|
||||
}
|
||||
|
||||
PathCompute->compute_peratom();
|
||||
double **path = PathCompute->array_atom;
|
||||
|
||||
double xum=0.;
|
||||
|
||||
// proj 0,1,2 = f.n, v.n, h.n
|
||||
// proj 3,4,5 = psi, f.n**2, f*(1-psi)
|
||||
// c_v 0,1,2 = fxcom, fycom, fzcom etc
|
||||
for(int i = 0; i < 10; i++) {
|
||||
c_v[i] = 0.;
|
||||
c_v_all[i] = 0.;
|
||||
}
|
||||
for(int i = 0; i < 6; i++) {
|
||||
proj[i] = 0.;
|
||||
proj_all[i] = 0.;
|
||||
}
|
||||
|
||||
double deviation[3] = {0.,0.,0.};
|
||||
|
||||
double fn;
|
||||
|
||||
force_flag=0;
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
h[i][0] = random->uniform() - 0.5;
|
||||
h[i][1] = random->uniform() - 0.5;
|
||||
h[i][2] = random->uniform() - 0.5;
|
||||
|
||||
proj[0] += f[i][0] * path[i][3]; // f.n
|
||||
proj[0] += f[i][1] * path[i][4]; // f.n
|
||||
proj[0] += f[i][2] * path[i][5]; // f.n
|
||||
|
||||
proj[1] += v[i][0] * path[i][3]; // v.n
|
||||
proj[1] += v[i][1] * path[i][4]; // v.n
|
||||
proj[1] += v[i][2] * path[i][5]; // v.n
|
||||
|
||||
proj[2] += h[i][0] * path[i][3]; // h.n
|
||||
proj[2] += h[i][1] * path[i][4]; // h.n
|
||||
proj[2] += h[i][2] * path[i][5]; // h.n
|
||||
|
||||
deviation[0] = x[i][0]-path[i][0]; // x-path
|
||||
deviation[1] = x[i][1]-path[i][1]; // x-path
|
||||
deviation[2] = x[i][2]-path[i][2]; // x-path
|
||||
domain->minimum_image(deviation);
|
||||
|
||||
proj[3] += path[i][6]*deviation[0]; // (x-path).dn/nn = psi
|
||||
proj[3] += path[i][7]*deviation[1]; // (x-path).dn/nn = psi
|
||||
proj[3] += path[i][8]*deviation[2]; // (x-path).dn/nn = psi
|
||||
|
||||
proj[4] += path[i][3]*deviation[0]; // (x-path).n
|
||||
proj[4] += path[i][4]*deviation[1]; // (x-path).n
|
||||
proj[4] += path[i][5]*deviation[2]; // (x-path).n
|
||||
|
||||
proj[5] += f[i][3]*deviation[0]; // (x-path).f
|
||||
proj[5] += f[i][4]*deviation[1]; // (x-path).f
|
||||
proj[5] += f[i][5]*deviation[2]; // (x-path).f
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
if(com_flag == 0){
|
||||
c_v[9] += 1.0;
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
c_v[0] += f[i][0];
|
||||
c_v[1] += f[i][1];
|
||||
c_v[2] += f[i][2];
|
||||
|
||||
c_v[3] += v[i][0];
|
||||
c_v[4] += v[i][1];
|
||||
c_v[5] += v[i][2];
|
||||
|
||||
c_v[6] += h[i][0];
|
||||
c_v[7] += h[i][1];
|
||||
c_v[8] += h[i][2];
|
||||
|
||||
c_v[9] += 1.0;
|
||||
}
|
||||
}
|
||||
MPI_Allreduce(proj,proj_all,6,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(c_v,c_v_all,10,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// results - f.n*(1-psi), (f.n)^2*(1-psi)^2, 1-psi, dX.n
|
||||
results_all[0] = proj_all[0] * (1.-proj_all[3]);
|
||||
results_all[1] = results_all[0] * results_all[0];
|
||||
results_all[2] = 1.-proj_all[3];
|
||||
results_all[3] = fabs(proj_all[4]);
|
||||
results_all[4] = proj_all[5]; // dX.f
|
||||
force_flag = 1;
|
||||
|
||||
for (int i = 0; i < nlocal; i++){
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
f[i][0] -= proj_all[0] * path[i][3] + c_v_all[0]/c_v_all[9];
|
||||
f[i][1] -= proj_all[0] * path[i][4] + c_v_all[1]/c_v_all[9];
|
||||
f[i][2] -= proj_all[0] * path[i][5] + c_v_all[2]/c_v_all[9];
|
||||
|
||||
v[i][0] -= proj_all[1] * path[i][3] + c_v_all[3]/c_v_all[9];
|
||||
v[i][1] -= proj_all[1] * path[i][4] + c_v_all[4]/c_v_all[9];
|
||||
v[i][2] -= proj_all[1] * path[i][5] + c_v_all[5]/c_v_all[9];
|
||||
|
||||
h[i][0] -= proj_all[2] * path[i][3] + c_v_all[6]/c_v_all[9];
|
||||
h[i][1] -= proj_all[2] * path[i][4] + c_v_all[7]/c_v_all[9];
|
||||
h[i][2] -= proj_all[2] * path[i][5] + c_v_all[8]/c_v_all[9];
|
||||
}
|
||||
}
|
||||
|
||||
if (od_flag == 0) {
|
||||
for (int i = 0; i < nlocal; i++){
|
||||
if (mask[i] & groupbit) {
|
||||
if(rmass) mass_f = sqrt(rmass[i]);
|
||||
else mass_f = sqrt(mass[type[i]]);
|
||||
|
||||
f[i][0] += -gamma * mass_f * mass_f * v[i][0];
|
||||
f[i][1] += -gamma * mass_f * mass_f * v[i][1];
|
||||
f[i][2] += -gamma * mass_f * mass_f * v[i][2];
|
||||
|
||||
f[i][0] += sqrtD * mass_f * h[i][0];
|
||||
f[i][1] += sqrtD * mass_f * h[i][1];
|
||||
f[i][2] += sqrtD * mass_f * h[i][2];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++){
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
if(rmass) mass_f = sqrt(rmass[i]);
|
||||
else mass_f = sqrt(mass[type[i]]);
|
||||
|
||||
f[i][0] += sqrtD * h[i][0] * mass_f;
|
||||
f[i][1] += sqrtD * h[i][1] * mass_f;
|
||||
f[i][2] += sqrtD * h[i][2] * mass_f;
|
||||
|
||||
f[i][0] /= gamma * mass_f * mass_f;
|
||||
f[i][1] /= gamma * mass_f * mass_f;
|
||||
f[i][2] /= gamma * mass_f * mass_f;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
void FixPAFI::post_force_respa(int vflag, int ilevel, int iloop)
|
||||
{
|
||||
// set force to desired value on requested level, 0.0 on other levels
|
||||
|
||||
if (ilevel == ilevel_respa) post_force(vflag);
|
||||
else {
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
for (int k = 0; k < 3; k++) f[i][k] = 0.0;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
void FixPAFI::min_post_force(int vflag)
|
||||
{
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double *rmass = atom->rmass;
|
||||
double *mass = atom->mass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
PathCompute->compute_peratom();
|
||||
double **path = PathCompute->array_atom;
|
||||
|
||||
double xum=0.;
|
||||
|
||||
// proj 0,1,2 = f.n, v.n, h.n
|
||||
// proj 3,4,5 = psi, f.n**2, f*(1-psi)
|
||||
// c_v 0,1,2 = fxcom, fycom, fzcom etc
|
||||
for(int i = 0; i < 10; i++) {
|
||||
c_v[i] = 0.;
|
||||
c_v_all[i] = 0.;
|
||||
}
|
||||
for(int i = 0; i < 6; i++) {
|
||||
proj[i] = 0.;
|
||||
proj_all[i] = 0.;
|
||||
}
|
||||
|
||||
double deviation[3] = {0.,0.,0.};
|
||||
|
||||
double fn;
|
||||
|
||||
force_flag=0;
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
proj[0] += f[i][0] * path[i][3]; // f.n
|
||||
proj[0] += f[i][1] * path[i][4]; // f.n
|
||||
proj[0] += f[i][2] * path[i][5]; // f.n
|
||||
|
||||
proj[1] += v[i][0] * path[i][3]; // v.n
|
||||
proj[1] += v[i][1] * path[i][4]; // v.n
|
||||
proj[1] += v[i][2] * path[i][5]; // v.n
|
||||
|
||||
proj[2] += h[i][0] * path[i][3]; // h.n
|
||||
proj[2] += h[i][1] * path[i][4]; // h.n
|
||||
proj[2] += h[i][2] * path[i][5]; // h.n
|
||||
|
||||
deviation[0] = x[i][0]-path[i][0]; // x-path
|
||||
deviation[1] = x[i][1]-path[i][1]; // x-path
|
||||
deviation[2] = x[i][2]-path[i][2]; // x-path
|
||||
domain->minimum_image(deviation);
|
||||
|
||||
proj[3] += path[i][6]*deviation[0]; // (x-path).dn/nn = psi
|
||||
proj[3] += path[i][7]*deviation[1]; // (x-path).dn/nn = psi
|
||||
proj[3] += path[i][8]*deviation[2]; // (x-path).dn/nn = psi
|
||||
|
||||
proj[4] += path[i][3]*deviation[0]; // (x-path).n
|
||||
proj[4] += path[i][4]*deviation[1]; // (x-path).n
|
||||
proj[4] += path[i][5]*deviation[2]; // (x-path).n
|
||||
|
||||
proj[5] += f[i][3]*deviation[0]; // (x-path).f
|
||||
proj[5] += f[i][4]*deviation[1]; // (x-path).f
|
||||
proj[5] += f[i][5]*deviation[2]; // (x-path).f
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
if(com_flag == 0){
|
||||
c_v[9] += 1.0;
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
c_v[0] += f[i][0];
|
||||
c_v[1] += f[i][1];
|
||||
c_v[2] += f[i][2];
|
||||
|
||||
c_v[3] += v[i][0];
|
||||
c_v[4] += v[i][1];
|
||||
c_v[5] += v[i][2];
|
||||
|
||||
c_v[6] += h[i][0];
|
||||
c_v[7] += h[i][1];
|
||||
c_v[8] += h[i][2];
|
||||
|
||||
c_v[9] += 1.0;
|
||||
}
|
||||
}
|
||||
MPI_Allreduce(proj,proj_all,6,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(c_v,c_v_all,10,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
results_all[0] = proj_all[0] * (1.-proj_all[3]); // f.n * psi
|
||||
results_all[1] = results_all[0] * results_all[0]; // (f.n * psi)^2
|
||||
results_all[2] = 1.-proj_all[3]; // psi
|
||||
results_all[3] = fabs(proj_all[4]); // dX.n
|
||||
results_all[4] = proj_all[5]; // dX.f
|
||||
|
||||
MPI_Bcast(results_all,5,MPI_DOUBLE,0,world);
|
||||
force_flag = 1;
|
||||
|
||||
for (int i = 0; i < nlocal; i++){
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
f[i][0] -= proj_all[0] * path[i][3] + c_v_all[0]/c_v_all[9];
|
||||
f[i][1] -= proj_all[0] * path[i][4] + c_v_all[1]/c_v_all[9];
|
||||
f[i][2] -= proj_all[0] * path[i][5] + c_v_all[2]/c_v_all[9];
|
||||
|
||||
v[i][0] -= proj_all[1] * path[i][3] + c_v_all[3]/c_v_all[9];
|
||||
v[i][1] -= proj_all[1] * path[i][4] + c_v_all[4]/c_v_all[9];
|
||||
v[i][2] -= proj_all[1] * path[i][5] + c_v_all[5]/c_v_all[9];
|
||||
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
double FixPAFI::compute_vector(int n)
|
||||
{
|
||||
return results_all[n];
|
||||
};
|
||||
|
||||
|
||||
|
||||
void FixPAFI::initial_integrate(int vflag)
|
||||
{
|
||||
double dtfm;
|
||||
|
||||
// update v and x of atoms in group
|
||||
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
|
||||
double *rmass = atom->rmass;
|
||||
double *mass = atom->mass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
PathCompute->compute_peratom();
|
||||
double **path = PathCompute->array_atom;
|
||||
|
||||
for(int i = 0; i < 10; i++) {
|
||||
c_v[i] = 0.;
|
||||
c_v_all[i] = 0.;
|
||||
}
|
||||
for(int i = 0; i < 6; i++) {
|
||||
proj[i] = 0.;
|
||||
proj_all[i] = 0.;
|
||||
}
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
proj[0] += f[i][0] * path[i][3]; // f.n
|
||||
proj[0] += f[i][1] * path[i][4]; // f.n
|
||||
proj[0] += f[i][2] * path[i][5]; // f.n
|
||||
|
||||
proj[1] += v[i][0] * path[i][3]; // v.n
|
||||
proj[1] += v[i][1] * path[i][4]; // v.n
|
||||
proj[1] += v[i][2] * path[i][5]; // v.n
|
||||
}
|
||||
}
|
||||
if(com_flag == 0){
|
||||
c_v[9] += 1.0;
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
c_v[0] += v[i][0];
|
||||
c_v[1] += v[i][1];
|
||||
c_v[2] += v[i][2];
|
||||
|
||||
c_v[3] += f[i][0];
|
||||
c_v[4] += f[i][1];
|
||||
c_v[5] += f[i][2];
|
||||
|
||||
c_v[9] += 1.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
MPI_Allreduce(proj,proj_all,5,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(c_v,c_v_all,10,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
if (od_flag == 0){
|
||||
if (rmass) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf / rmass[i];
|
||||
v[i][0] += dtfm * (f[i][0]-path[i][3]*proj_all[0] - c_v_all[3]/c_v_all[9]);
|
||||
v[i][1] += dtfm * (f[i][1]-path[i][4]*proj_all[0] - c_v_all[4]/c_v_all[9]);
|
||||
v[i][2] += dtfm * (f[i][2]-path[i][5]*proj_all[0] - c_v_all[5]/c_v_all[9]);
|
||||
x[i][0] += dtv * (v[i][0]-path[i][3]*proj_all[1] - c_v_all[0]/c_v_all[9]);
|
||||
x[i][1] += dtv * (v[i][1]-path[i][4]*proj_all[1] - c_v_all[1]/c_v_all[9]);
|
||||
x[i][2] += dtv * (v[i][2]-path[i][5]*proj_all[1] - c_v_all[2]/c_v_all[9]);
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf / mass[type[i]];
|
||||
v[i][0] += dtfm * (f[i][0]-path[i][3]*proj_all[0] - c_v_all[3]/c_v_all[9]);
|
||||
v[i][1] += dtfm * (f[i][1]-path[i][4]*proj_all[0] - c_v_all[4]/c_v_all[9]);
|
||||
v[i][2] += dtfm * (f[i][2]-path[i][5]*proj_all[0] - c_v_all[5]/c_v_all[9]);
|
||||
x[i][0] += dtv * (v[i][0]-path[i][3]*proj_all[1] - c_v_all[0]/c_v_all[9]);
|
||||
x[i][1] += dtv * (v[i][1]-path[i][4]*proj_all[1] - c_v_all[1]/c_v_all[9]);
|
||||
x[i][2] += dtv * (v[i][2]-path[i][5]*proj_all[1] - c_v_all[2]/c_v_all[9]);
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if (rmass) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf / rmass[i];
|
||||
v[i][0] = 0.;
|
||||
v[i][1] = 0.;
|
||||
v[i][2] = 0.;
|
||||
x[i][0] += dtv * (f[i][0]-path[i][3]*proj_all[0] - c_v_all[3]/c_v_all[9]);
|
||||
x[i][1] += dtv * (f[i][1]-path[i][4]*proj_all[0] - c_v_all[4]/c_v_all[9]);
|
||||
x[i][2] += dtv * (f[i][2]-path[i][5]*proj_all[0] - c_v_all[5]/c_v_all[9]);
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf / mass[type[i]];
|
||||
v[i][0] = 0.;
|
||||
v[i][1] = 0.;
|
||||
v[i][2] = 0.;
|
||||
x[i][0] += dtv * (f[i][0]-path[i][3]*proj_all[0] - c_v_all[3]/c_v_all[9]);
|
||||
x[i][1] += dtv * (f[i][1]-path[i][4]*proj_all[0] - c_v_all[4]/c_v_all[9]);
|
||||
x[i][2] += dtv * (f[i][2]-path[i][5]*proj_all[0] - c_v_all[5]/c_v_all[9]);
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixPAFI::final_integrate()
|
||||
{
|
||||
double dtfm;
|
||||
|
||||
// update v of atoms in group
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double *rmass = atom->rmass;
|
||||
double *mass = atom->mass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
PathCompute->compute_peratom();
|
||||
double **path = PathCompute->array_atom;
|
||||
|
||||
for(int i = 0; i < 10; i++) {
|
||||
c_v[i] = 0.;
|
||||
c_v_all[i] = 0.;
|
||||
}
|
||||
for(int i = 0; i < 6; i++) {
|
||||
proj[i] = 0.;
|
||||
proj_all[i] = 0.;
|
||||
}
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
proj[0] += f[i][0] * path[i][3]; // f.n
|
||||
proj[0] += f[i][1] * path[i][4]; // f.n
|
||||
proj[0] += f[i][2] * path[i][5]; // f.n
|
||||
}
|
||||
if(com_flag == 0){
|
||||
c_v[9] += 1.0;
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
c_v[3] += f[i][0];
|
||||
c_v[4] += f[i][1];
|
||||
c_v[5] += f[i][2];
|
||||
c_v[9] += 1.0;
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(proj,proj_all,5,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(c_v,c_v_all,10,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
if (od_flag == 0){
|
||||
if (rmass) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf / rmass[i];
|
||||
v[i][0] += dtfm * (f[i][0]-path[i][3]*proj_all[0] - c_v_all[3]/c_v_all[9]);
|
||||
v[i][1] += dtfm * (f[i][1]-path[i][4]*proj_all[0] - c_v_all[4]/c_v_all[9]);
|
||||
v[i][2] += dtfm * (f[i][2]-path[i][5]*proj_all[0] - c_v_all[5]/c_v_all[9]);
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf / mass[type[i]];
|
||||
v[i][0] += dtfm * (f[i][0]-path[i][3]*proj_all[0] - c_v_all[3]/c_v_all[9]);
|
||||
v[i][1] += dtfm * (f[i][1]-path[i][4]*proj_all[0] - c_v_all[4]/c_v_all[9]);
|
||||
v[i][2] += dtfm * (f[i][2]-path[i][5]*proj_all[0] - c_v_all[5]/c_v_all[9]);
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
v[i][0] = 0.;
|
||||
v[i][1] = 0.;
|
||||
v[i][2] = 0.;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixPAFI::initial_integrate_respa(int vflag, int ilevel, int iloop)
|
||||
{
|
||||
dtv = step_respa[ilevel];
|
||||
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
|
||||
|
||||
// innermost level - NVE update of v and x
|
||||
// all other levels - NVE update of v
|
||||
|
||||
if (ilevel == 0) initial_integrate(vflag);
|
||||
else final_integrate();
|
||||
};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixPAFI::final_integrate_respa(int ilevel, int iloop)
|
||||
{
|
||||
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
|
||||
final_integrate();
|
||||
};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixPAFI::reset_dt()
|
||||
{
|
||||
dtv = update->dt;
|
||||
dtf = 0.5 * update->dt * force->ftm2v;
|
||||
};
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of local atom-based array
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixPAFI::memory_usage()
|
||||
{
|
||||
double bytes = 0.0;
|
||||
bytes = maxatom*3 * sizeof(double);
|
||||
return bytes;
|
||||
};
|
||||
|
|
@ -0,0 +1,99 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(pafi,FixPAFI)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_PAFI_H
|
||||
#define LMP_FIX_PAFI_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixPAFI : public Fix {
|
||||
public:
|
||||
FixPAFI(class LAMMPS *, int, char **);
|
||||
virtual ~FixPAFI();
|
||||
int setmask();
|
||||
virtual void init();
|
||||
|
||||
void setup(int);
|
||||
void min_setup(int);
|
||||
virtual void post_force(int);
|
||||
|
||||
void post_force_respa(int, int, int);
|
||||
void min_post_force(int);
|
||||
double compute_vector(int);
|
||||
// nve
|
||||
virtual void initial_integrate(int);
|
||||
virtual void final_integrate();
|
||||
virtual void initial_integrate_respa(int, int, int);
|
||||
virtual void final_integrate_respa(int, int);
|
||||
virtual void reset_dt();
|
||||
|
||||
double memory_usage();
|
||||
|
||||
protected:
|
||||
int varflag,icompute;
|
||||
char *computename;
|
||||
class Compute *PathCompute;
|
||||
double proj[6], proj_all[6]; // f,v,h, psi
|
||||
double results[5], results_all[5]; // f.n, (f.n)**2, psi, dx.n
|
||||
double c_v[10],c_v_all[10];
|
||||
double temperature,gamma,sqrtD,t_period,local_norm,mass_f;
|
||||
int force_flag,od_flag,com_flag;
|
||||
int nlevels_respa,ilevel_respa;
|
||||
int maxatom;
|
||||
class RanMars *random;
|
||||
int seed;
|
||||
double **h;
|
||||
// nve
|
||||
double dtv,dtf;
|
||||
double *step_respa;
|
||||
int mass_require;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Compute ID for fix pafi does not exist
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Compute for fix pafi does not calculate a local array
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Compute for fix pafi has < 3*DIM fields per atom
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix pafi requires a damped dynamics minimizer
|
||||
|
||||
Use a different minimization style.
|
||||
|
||||
*/
|
||||
1128
src/library.cpp
1128
src/library.cpp
File diff suppressed because it is too large
Load Diff
|
|
@ -132,6 +132,14 @@ int lammps_set_variable(void *, char *, char *);
|
|||
* Library functions for scatter/gather operations of data
|
||||
* ---------------------------------------------------------------------- */
|
||||
|
||||
|
||||
void lammps_gather(void *, char *, int, int, void *);
|
||||
void lammps_gather_concat(void *, char *, int, int, void *);
|
||||
void lammps_gather_subset(void *, char *, int, int, int, int *, void *);
|
||||
void lammps_scatter(void *, char *, int, int, void *);
|
||||
void lammps_scatter_subset(void *, char *, int, int, int, int *, void *);
|
||||
|
||||
|
||||
void lammps_gather_atoms(void *, char *, int, int, void *);
|
||||
void lammps_gather_atoms_concat(void *, char *, int, int, void *);
|
||||
void lammps_gather_atoms_subset(void *, char *, int, int, int, int *, void *);
|
||||
|
|
|
|||
Loading…
Reference in New Issue