Skip to content
Snippets Groups Projects
Commit e53583d9 authored by sjplimp's avatar sjplimp Committed by GitHub
Browse files

Merge pull request #590 from lammps/fortran-dftb 

update of Fortran-DFTB interface to be compatible with fix msst
parents 3f833968 551001f1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 344 additions and 43 deletions
doc/src/fix_msst.txt
+ 27
− 19
View file @ e53583d9
......@@ -25,7 +25,7 @@ keyword = {q} or {mu} or {p0} or {v0} or {e0} or {tscale} or {beta} or {dftb} :l
{e0} value = initial total energy (energy units)
{tscale} value = reduction in initial temperature (unitless fraction between 0.0 and 1.0)
{dftb} value = {yes} or {no} for whether using MSST in conjunction with DFTB+
{beta} value = scale factor on energy contribution of DFTB+ :pre
{beta} value = scale factor for improved energy conservation :pre
:ule
[Examples:]
......@@ -72,6 +72,14 @@ be calculated on the first step, after the energy specified by
{tscale} is removed. The value of {e0} is not used in the dynamical
equations, but is used in calculating the deviation from the Hugoniot.
The keyword {beta} is a scaling term that can be added to the MSST
ionic equations of motion to account for drift in the conserved
quantity during long timescale simulations, similar to a Berendson
thermostat. See "(Reed)"_#Reed and "(Goldman)"_#Goldman for more
details. The value of {beta} must be between 0.0 and 1.0 inclusive.
A value of 0.0 means no contribution, a value of 1.0 means a full
contribution.
Values of shockvel less than a critical value determined by the
material response will not have compressive solutions. This will be
reflected in lack of significant change of the volume in the MSST.
......@@ -95,23 +103,15 @@ or "_MSST_pe". The group for the new computes is "all".
:line
The {dftb} and {beta} keywords are to allow this fix to be used when
LAMMPS is being driven by DFTB+, a density-functional tight-binding
code.
If the keyword {dftb} is used with a value of {yes}, then the MSST
equations are altered to account for an energy contribution compute by
DFTB+. In this case, you must define a "fix
external"_fix_external.html command in your input script, which is
used to callback to DFTB+ during the LAMMPS timestepping. DFTB+ will
communicate its info to LAMMPS via that fix.
The keyword {beta} is a scale factor on the DFTB+ energy contribution.
The value of {beta} must be between 0.0 and 1.0 inclusive. A value of
0.0 means no contribution, a value of 1.0 means a full contribution.
(July 2017) More information about these keywords and the use of
LAMMPS with DFTB+ will be added to the LAMMMPS documention soon.
The {dftb} keyword is to allow this fix to be used when LAMMPS is
being driven by DFTB+, a density-functional tight-binding code. If the
keyword {dftb} is used with a value of {yes}, then the MSST equations
are altered to account for the electron entropy contribution to the
Hugonio relations and total energy. See "(Reed2)"_#Reed2 and
"(Goldman)"_#Goldman for details on this contribution. In this case,
you must define a "fix external"_fix_external.html command in your
input script, which is used to callback to DFTB+ during the LAMMPS
timestepping. DFTB+ will communicate its info to LAMMPS via that fix.
:line
......@@ -182,4 +182,12 @@ timestep.
:line
:link(Reed)
[(Reed)] Reed, Fried, and Joannopoulos, Phys. Rev. Lett., 90, 235503 (2003).
[(Reed)] Reed, Fried, and Joannopoulos, Phys. Rev. Lett., 90, 235503
(2003).
:link(Reed2)
[(Reed2)] Reed, J. Phys. Chem. C, 116, 2205 (2012).
:link(Goldman)
[(Goldman)] Goldman, Srinivasan, Hamel, Fried, Gaus, and Elstner,
J. Phys. Chem. C, 117, 7885 (2013).
examples/COUPLE/README
+ 2
− 2
View file @ e53583d9
......@@ -41,8 +41,8 @@ fortran a simple wrapper on the LAMMPS library API that
can be called from Fortran
fortran2 a more sophisticated wrapper on the LAMMPS library API that
can be called from Fortran
fortran3 wrapper written by Nir Goldman (LLNL), as an
fortran_dftb wrapper written by Nir Goldman (LLNL), as an
extension to fortran2, used for calling LAMMPS
from Fortran DFTB+ code
from Fortran DFTB+ tight-binding code
Each sub-directory has its own README with more details.
examples/COUPLE/fortran3/LAMMPS-wrapper2.cpp → examples/COUPLE/fortran_dftb/LAMMPS-wrapper2.cpp
+ 25
− 1
View file @ e53583d9
......@@ -47,11 +47,35 @@ void lammps_set_callback (void *ptr) {
return;
}
void lammps_set_external_vector_length (void *ptr, int n) {
class LAMMPS *lmp = (class LAMMPS *) ptr;
int ifix = lmp->modify->find_fix_by_style("external");
FixExternal *fix = (FixExternal *) lmp->modify->fix[ifix];
fix->set_vector_length(n);
return;
}
void lammps_set_external_vector (void *ptr, int n, double val) {
class LAMMPS *lmp = (class LAMMPS *) ptr;
int ifix = lmp->modify->find_fix_by_style("external");
FixExternal *fix = (FixExternal *) lmp->modify->fix[ifix];
fix->set_vector (n, val);
return;
}
void lammps_set_user_energy (void *ptr, double energy) {
class LAMMPS *lmp = (class LAMMPS *) ptr;
int ifix = lmp->modify->find_fix_by_style("external");
FixExternal *fix = (FixExternal *) lmp->modify->fix[ifix];
fix->set_energy(energy);
fix->set_energy_global(energy);
return;
}
void lammps_set_user_virial (void *ptr, double *virial) {
class LAMMPS *lmp = (class LAMMPS *) ptr;
int ifix = lmp->modify->find_fix_by_style("external");
FixExternal *fix = (FixExternal *) lmp->modify->fix[ifix];
fix->set_virial_global(virial);
return;
}
examples/COUPLE/fortran3/LAMMPS-wrapper2.h → examples/COUPLE/fortran_dftb/LAMMPS-wrapper2.h
+ 3
− 0
View file @ e53583d9
......@@ -26,6 +26,9 @@ extern "C" {
/* Prototypes for auxiliary functions */
void lammps_set_callback (void *);
void lammps_set_user_energy (void*, double);
void lammps_set_user_virial (void*, double*);
void lammps_set_external_vector_length (void*, int);
void lammps_set_external_vector (void*, int, double);
#ifdef __cplusplus
}
......
examples/COUPLE/fortran3/LAMMPS.F90 → examples/COUPLE/fortran_dftb/LAMMPS.F90
+ 28
− 2
View file @ e53583d9
......@@ -52,12 +52,16 @@ module LAMMPS
C_NULL_CHAR, C_loc, C_F_pointer, lammps_instance => C_ptr
implicit none
private
public :: lammps_set_user_virial
public :: lammps_set_external_vector_length
public :: lammps_set_external_vector
public :: lammps_set_user_energy
public :: lammps_open, lammps_open_no_mpi, lammps_close, lammps_file, &
lammps_command, lammps_free, lammps_extract_global, &
lammps_extract_atom, lammps_extract_compute, lammps_extract_fix, &
lammps_extract_variable, lammps_get_natoms, lammps_gather_atoms, &
lammps_scatter_atoms, lammps_set_callback, lammps_set_user_energy
public :: lammps_instance, C_ptr, C_double, C_int
lammps_set_callback
public :: lammps_scatter_atoms, lammps_instance, C_ptr, C_double, C_int
!! Functions supplemental to the prototypes in library.h. {{{1
!! The function definitions (in C++) are contained in LAMMPS-wrapper.cpp.
......@@ -218,6 +222,28 @@ module LAMMPS
real(C_double), value :: energy
end subroutine lammps_set_user_energy
subroutine lammps_set_user_virial (ptr, virial) &
bind (C, name='lammps_set_user_virial')
import :: C_ptr, C_double
type (C_ptr), value :: ptr
real(C_double) :: virial(6)
end subroutine lammps_set_user_virial
subroutine lammps_set_external_vector_length (ptr, n) &
bind (C, name='lammps_set_external_vector_length')
import :: C_ptr, C_double, C_int
type(C_ptr), value :: ptr
integer (C_int), value :: n
end subroutine lammps_set_external_vector_length
subroutine lammps_set_external_vector (ptr, n, val) &
bind (C, name='lammps_set_external_vector')
import :: C_ptr, C_int, C_double
type (C_ptr), value :: ptr
integer (C_int), value :: n
real(C_double), value :: val
end subroutine lammps_set_external_vector
subroutine lammps_actual_gather_atoms (ptr, name, type, count, data) &
bind (C, name='lammps_gather_atoms')
import :: C_ptr, C_int, C_char
......
examples/COUPLE/fortran3/README → examples/COUPLE/fortran_dftb/README
+ 8
− 4
View file @ e53583d9
......@@ -3,8 +3,9 @@ forces from a fortran code for a LAMMPS simulation. The reader should
refer to the README file in COUPLE/fortran2 before proceeding. Here,
the LAMMPS.F90 file has been modified slightly and additional files
named LAMMPS-wrapper2.h and LAMMPS-wrapper2.cpp have been included in
order to supply wrapper functions to set the LAMMPS callback function
and total energy.
order to supply wrapper functions to set the LAMMPS callback function,
total energy, virial, and electronic entropy contribution (needed for
MSST simulations with a quantum code).
In this example, the callback function is set to run the
semi-empirical quantum code DFTB+ in serial and then read in the total
......@@ -20,11 +21,14 @@ etc.
A few more important notes:
-The stress tensor from DFTB+ is passed in to LAMMPS via pointer.
-Calling the subroutine lammps_set_callback() is required in order to set
a pointer to the callback function in LAMMPS.
-The subroutine lammps_set_user_energy() passes in the potential energy
from DFTB+ to LAMMPS.
from DFTB+ to LAMMPS. Similarly, lammps_set_user_virial passes the stress tensor.
-The electronic entropy contribution is set via lammps_set_external_vector(). Their needs
to be a call to lammps_set_external_vector_length() before this value can be
passed to LAMMPS.
This example was created by Nir Goldman, whom you can contact with
questions:
......
examples/COUPLE/fortran_dftb/dftb_in.hsd 0 → 100644
+ 40
− 0
View file @ e53583d9
#sample DFTB+ script to run this test code
Geometry = GenFormat {
<<< "lammps.gen"
}
Driver = {
}
Hamiltonian = DFTB {
LAMMPS = Yes # keyword to print energy, forces, and stress tensor to file(results.out)
SCC = No
MaxAngularMomentum = {
C = "p"
}
Charge = 0.0
Eigensolver = Standard {}
Filling = Fermi {
Temperature [Kelvin] = 298.0
}
SlaterKosterFiles = Type2FileNames {
Prefix = "~/slako/mio-1-1/" # the user must define the location of the skf files
Separator = "-"
Suffix = ".skf"
LowerCaseTypeName = No
}
KPointsAndWeights = {
0.0000000000000 0.0000000000000 0.0000000000000 1.00000000000000
}
}
Options = {
CalculateForces = Yes
WriteDetailedOut = No
WriteBandOut = No
RandomSeed = 12345
}
ParserOptions = {
ParserVersion = 3
}
examples/COUPLE/fortran_dftb/dftb_pin.hsd 0 → 100644
+ 129
− 0
View file @ e53583d9
Geometry = GenFormat {
64 S
C
1 1 7.099007 7.117657 7.119139
2 1 0.858709 0.867233 0.882294
3 1 1.772527 1.811776 7.120239
4 1 2.702145 2.681271 0.901362
5 1 0.017539 1.794455 1.788454
6 1 0.885593 2.694118 2.707994
7 1 1.795055 7.120787 1.777896
8 1 2.642849 0.868278 2.670699
9 1 0.016060 0.017156 3.568644
10 1 0.891891 0.896406 4.439286
11 1 1.766086 1.764402 3.550134
12 1 2.677349 2.648926 4.427174
13 1 0.010133 1.771283 5.342173
14 1 0.858153 2.653565 6.241596
15 1 1.804087 0.020636 5.353268
16 1 2.689680 0.907188 6.224575
17 1 0.017845 3.577563 7.113016
18 1 0.910027 4.459286 0.910286
19 1 1.766394 5.376046 0.015526
20 1 2.683727 6.220728 0.898553
21 1 0.003357 5.363423 1.774139
22 1 0.856735 6.238324 2.660213
23 1 1.761079 3.549776 1.797054
24 1 2.667227 4.463441 2.646074
25 1 7.132499 3.551558 3.599764
26 1 0.920387 4.482191 4.479257
27 1 1.772194 5.337132 3.555569
28 1 2.675010 6.251629 4.483124
29 1 0.005702 5.371095 5.351147
30 1 0.880807 6.249819 6.264231
31 1 1.793177 3.592396 5.369939
32 1 2.653179 4.463595 6.274044
33 1 3.557243 7.118913 0.026006
34 1 4.458971 0.889331 0.904950
35 1 5.367903 1.759757 7.104941
36 1 6.271565 2.658454 0.890168
37 1 3.591915 1.768681 1.793880
38 1 4.435612 2.662184 2.676722
39 1 5.371040 0.000196 1.783464
40 1 6.226453 0.886640 2.653384
41 1 3.583339 0.005449 3.600177
42 1 4.453692 0.909417 4.459713
43 1 5.314554 1.805409 3.584215
44 1 6.210181 2.642660 4.486206
45 1 3.545704 1.802745 5.365369
46 1 4.476660 2.701226 6.220451
47 1 5.332820 0.029557 5.347965
48 1 6.215725 0.915081 6.230289
49 1 3.536446 3.551469 7.106600
50 1 4.451181 4.426439 0.900180
51 1 5.368735 5.377996 7.109524
52 1 6.230666 6.220985 0.862175
53 1 3.596626 5.372822 1.797613
54 1 4.485613 6.221252 2.699652
55 1 5.364421 3.549838 1.796281
56 1 6.261739 4.459046 2.648152
57 1 3.588752 3.581054 3.581755
58 1 4.462342 4.467270 4.478800
59 1 5.355202 5.318323 3.556531
60 1 6.268570 6.259831 4.465795
61 1 3.588636 5.354278 5.362327
62 1 4.475747 6.263866 6.227803
63 1 5.331158 3.554349 5.318368
64 1 6.254581 4.436344 6.209681
0.0 0.0 0.0
7.13400000000000 0 0
0 7.13400000000000 0
0 0 7.13400000000000
}
Driver = {}
Hamiltonian = DFTB {
LAMMPS = Yes
SCC = No
MaxAngularMomentum = {
C = "p"
}
Charge = 0.0
Eigensolver = Standard {}
Filling = Fermi {
Temperature [Kelvin] = 298.0
IndependentKFilling = No
}
SlaterKosterFiles = Type2FileNames {
Prefix = "~/slako/mio-1-1/"
Separator = "-"
Suffix = ".skf"
LowerCaseTypeName = No
}
KPointsAndWeights = {
0.0000000000000 0.0000000000000 0.0000000000000 1.00000000000000
}
PolynomialRepulsive = {}
OldRepulsiveSum = No
OrbitalResolvedSCC = No
OldSKInterpolation = No
NoErep = No
Dispersion = {}
ThirdOrder = No
ThirdOrderFull = No
}
Options = {
CalculateForces = Yes
WriteDetailedOut = No
WriteBandOut = No
RandomSeed = 12345
MullikenAnalysis = No
WriteEigenvectors = No
WriteAutotestTag = No
WriteDetailedXML = No
WriteResultsTag = No
AtomResolvedEnergies = No
WriteHS = No
WriteRealHS = No
MinimiseMemoryUsage = No
ShowFoldedCoords = No
}
ParserOptions = {
ParserVersion = 3
WriteHSDInput = Yes
WriteXMLInput = No
StopAfterParsing = No
IgnoreUnprocessedNodes = No
}
Analysis = {
ProjectStates = {}
}
examples/COUPLE/fortran_dftb/log.simple 0 → 100644
+ 71
− 0
View file @ e53583d9
LAMMPS (6 Jul 2017)
units real
atom_style charge
atom_modify map array
atom_modify sort 0 0.0
read_data data.diamond
triclinic box = (0 0 0) to (7.134 7.134 7.134) with tilt (0 0 0)
1 by 1 by 1 MPI processor grid
reading atoms ...
64 atoms
reading velocities ...
64 velocities
neighbor 1.0 bin
neigh_modify delay 0 every 5 check no
fix 1 all nve
fix 2 all external pf/callback 1 1
fix_modify 2 energy yes
thermo_style custom step temp etotal ke pe lx ly lz pxx pyy pzz press
thermo 1
timestep 0.5
run 10
Neighbor list info ...
update every 5 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 0
ghost atom cutoff = 0
binsize = 7.134, bins = 1 1 1
0 neighbor lists, perpetual/occasional/extra = 0 0 0
Per MPI rank memory allocation (min/avg/max) = 2.3 | 2.3 | 2.3 Mbytes
Step Temp TotEng KinEng PotEng Lx Ly Lz Pxx Pyy Pzz Press
0 298.24835 -69593.587 56.008365 -69649.595 7.134 7.134 7.134 -19980.19 -21024.038 -21097.458 -20700.562
1 295.24358 -69593.585 55.444098 -69649.029 7.134 7.134 7.134 -19778.833 -20799.657 -20854.156 -20477.549
2 286.37211 -69593.58 53.778115 -69647.358 7.134 7.134 7.134 -19227.52 -20177.28 -20176.12 -19860.306
3 272.062 -69593.572 51.090804 -69644.663 7.134 7.134 7.134 -18360.869 -19189.684 -19100.021 -18883.525
4 253.01834 -69593.561 47.514575 -69641.075 7.134 7.134 7.134 -17198.143 -17855.03 -17652.036 -17568.403
5 230.19242 -69593.547 43.228073 -69636.775 7.134 7.134 7.134 -15750.247 -16183.764 -15854.145 -15929.386
6 204.71787 -69593.533 38.44418 -69631.977 7.134 7.134 7.134 -14083.498 -14247.434 -13789.835 -14040.256
7 177.82397 -69593.518 33.393748 -69626.911 7.134 7.134 7.134 -12340.963 -12202.878 -11623.171 -12055.671
8 150.76736 -69593.503 28.312758 -69621.816 7.134 7.134 7.134 -10637.824 -10180.827 -9495.0496 -10104.567
9 124.7737 -69593.49 23.431383 -69616.921 7.134 7.134 7.134 -9113.3842 -8339.0492 -7572.8076 -8341.747
10 100.98183 -69593.478 18.963481 -69612.442 7.134 7.134 7.134 -7833.9349 -6756.9749 -5945.8968 -6845.6022
Loop time of 2.20497 on 1 procs for 10 steps with 64 atoms
Performance: 0.196 ns/day, 122.499 hours/ns, 4.535 timesteps/s
0.2% 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 | 1.4305e-06 | 1.4305e-06 | 1.4305e-06 | 0.0 | 0.00
Comm | 4.22e-05 | 4.22e-05 | 4.22e-05 | 0.0 | 0.00
Output | 0.00067687 | 0.00067687 | 0.00067687 | 0.0 | 0.03
Modify | 2.2042 | 2.2042 | 2.2042 | 0.0 | 99.96
Other | | 6.533e-05 | | | 0.00
Nlocal: 64 ave 64 max 64 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 0
Ave neighs/atom = 0
Neighbor list builds = 2
Dangerous builds not checked
Total wall time: 0:00:02
examples/COUPLE/fortran3/makefile → examples/COUPLE/fortran_dftb/makefile
+ 1
− 1
View file @ e53583d9
......@@ -21,7 +21,7 @@ liblammps_fortran.so : LAMMPS.o LAMMPS-wrapper.o LAMMPS-wrapper2.o
$(FC) $(FFLAGS) -shared -o $@ $^
simpleF.x: simple.o LAMMPS.o LAMMPS-wrapper.o LAMMPS-wrapper2.o
$(FC) $(FFLAGS) simple.o -o simpleF.x liblammps_fortran.a $(LAMMPS_SRC)/liblammps_mvapich.a -lstdc++ /usr/local/tools/fftw/lib/libfftw.a
$(FC) $(FFLAGS) simple.o -o simpleF.x liblammps_fortran.a $(LAMMPS_SRC)/liblammps_mvapich.a -lstdc++ /usr/lib64/libfftw3.a
liblammps_fortran.a : LAMMPS.o LAMMPS-wrapper.o LAMMPS-wrapper2.o
$(AR) rs $@ $^
......
examples/COUPLE/fortran3/simple.f90 → examples/COUPLE/fortran_dftb/simple.f90
+ 10
− 14
View file @ e53583d9
......@@ -13,7 +13,7 @@
type(c_ptr) :: c_pos, c_fext, c_ids
double precision, pointer :: fext(:,:), pos(:,:)
integer, intent(in) :: ids(nlocal)
real (C_double), dimension(:), pointer :: virial => NULL()
real(C_double) :: virial(6)
real (C_double) :: etot
real(C_double), pointer :: ts_lmp
double precision :: stress(3,3), ts_dftb
......@@ -61,26 +61,21 @@
read(10,*)stress(i,:)
enddo
stress (:,:) = stress(:,:)*autoatm
virial(1) = stress(1,1)/(nktv2p/volume)
virial(2) = stress(2,2)/(nktv2p/volume)
virial(3) = stress(3,3)/(nktv2p/volume)
virial(4) = stress(1,2)/(nktv2p/volume)
virial(5) = stress(1,3)/(nktv2p/volume)
virial(6) = stress(2,3)/(nktv2p/volume)
etot = etot*econv
call lammps_extract_global(ts_lmp, lmp, 'TS_dftb')
ts_lmp = ts_dftb
call lammps_set_external_vector(lmp,1,ts_dftb*econv)
do i = 1, nlocal
read(10,*)fext(:,ids(i))
fext(:,ids(i)) = fext(:,ids(i))*fconv
enddo
close(10)
call lammps_set_user_energy (lmp, etot)
call lammps_extract_atom (virial, lmp, 'virial')
if (.not. associated(virial)) then
print*,'virial pointer not associated.'
STOP
endif
virial(1) = stress(1,1)/(nktv2p/volume)
virial(2) = stress(2,2)/(nktv2p/volume)
virial(3) = stress(3,3)/(nktv2p/volume)
virial(4) = stress(1,2)/(nktv2p/volume)
virial(5) = stress(1,3)/(nktv2p/volume)
virial(6) = stress(2,3)/(nktv2p/volume)
call lammps_set_user_virial (lmp, virial)
end subroutine
end module callback
......@@ -103,6 +98,7 @@ program simple_fortran_callback
call lammps_open_no_mpi ('lmp -log log.simple', lmp)
call lammps_file (lmp, 'in.simple')
call lammps_set_callback(lmp)
call lammps_set_external_vector_length(lmp,2)
call lammps_command (lmp, 'run 10')
call lammps_close (lmp)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment