2 |
|
!! Corresponds to the force field defined in lj_FF.cpp |
3 |
|
!! @author Charles F. Vardeman II |
4 |
|
!! @author Matthew Meineke |
5 |
< |
!! @version $Id: calc_LJ_FF.F90,v 1.11 2003-03-17 21:07:50 gezelter Exp $, $Date: 2003-03-17 21:07:50 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $ |
5 |
> |
!! @version $Id: calc_LJ_FF.F90,v 1.12 2003-03-18 16:46:47 gezelter Exp $, $Date: 2003-03-18 16:46:47 $, $Name: not supported by cvs2svn $, $Revision: 1.12 $ |
6 |
|
|
7 |
|
module lj |
8 |
– |
use simulation |
8 |
|
use definitions |
9 |
|
use atype_module |
10 |
|
use vector_class |
11 |
|
#ifdef IS_MPI |
12 |
|
use mpiSimulation |
13 |
|
#endif |
14 |
+ |
use force_globals |
15 |
+ |
|
16 |
|
implicit none |
17 |
|
PRIVATE |
18 |
|
|
20 |
|
#include "fForceField.h" |
21 |
|
|
22 |
|
integer, save :: LJ_Mixing_Policy |
23 |
+ |
integer, save :: LJ_rcut |
24 |
|
|
25 |
|
!! Logical has lj force field module been initialized? |
26 |
|
|
28 |
|
|
29 |
|
!! Public methods and data |
30 |
|
public :: init_LJ_FF |
31 |
+ |
public :: LJ_new_rcut |
32 |
|
public :: do_lj_pair |
33 |
|
|
34 |
|
type :: lj_mixed_params |
50 |
|
|
51 |
|
contains |
52 |
|
|
53 |
< |
subroutine init_LJ_FF(mix_Policy, status) |
53 |
> |
subroutine init_LJ_FF(mix_Policy, rcut, status) |
54 |
|
integer, intent(in) :: mix_Policy |
55 |
|
integer, intent(out) :: status |
56 |
|
integer :: myStatus |
57 |
+ |
real(kind=dp) :: rcut |
58 |
|
|
59 |
|
if (mix_Policy == LB_MIXING_RULE) then |
60 |
|
LJ_Mixing_Policy = LB_MIXING_RULE |
68 |
|
endif |
69 |
|
endif |
70 |
|
|
71 |
+ |
LJ_rcut = rcut |
72 |
+ |
|
73 |
|
status = 0 |
74 |
|
call createMixingList(myStatus) |
75 |
|
if (myStatus /= 0) then |
81 |
|
|
82 |
|
end subroutine init_LJ_FF |
83 |
|
|
84 |
+ |
subroutine LJ_new_rcut(rcut, status) |
85 |
+ |
integer :: status, myStatus |
86 |
+ |
real(kind=dp) :: rcut |
87 |
+ |
|
88 |
+ |
LJ_rcut = rcut |
89 |
+ |
status = 0 |
90 |
+ |
call createMixingList(myStatus) |
91 |
+ |
if (myStatus /= 0) then |
92 |
+ |
status = -1 |
93 |
+ |
return |
94 |
+ |
end if |
95 |
+ |
|
96 |
+ |
return |
97 |
+ |
end subroutine LJ_new_rcut |
98 |
+ |
|
99 |
|
subroutine createMixingList(status) |
100 |
|
integer :: nAtypes |
101 |
|
integer :: status |
103 |
|
integer :: j |
104 |
|
real ( kind = dp ) :: mySigma_i,mySigma_j |
105 |
|
real ( kind = dp ) :: myEpsilon_i,myEpsilon_j |
106 |
< |
real ( kind = dp ) :: rcut, rcut6 |
106 |
> |
real ( kind = dp ) :: rcut6 |
107 |
|
status = 0 |
108 |
|
|
109 |
|
nAtypes = getSize(atypes) |
118 |
|
status = -1 |
119 |
|
return |
120 |
|
end if |
121 |
< |
call getRcut(rcut, rc6=rcut6) |
121 |
> |
|
122 |
> |
rcut6 = LJ_rcut**6 |
123 |
> |
|
124 |
|
do i = 1, nAtypes |
125 |
|
|
126 |
< |
call getElementProperty(atypes,i,"lj_epsilon",myEpsilon_i) |
127 |
< |
call getElementProperty(atypes,i,"lj_sigma", mySigma_i) |
126 |
> |
call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i) |
127 |
> |
call getElementProperty(atypes, i, "lj_sigma", mySigma_i) |
128 |
|
! do self mixing rule |
129 |
|
ljMixed(i,i)%sigma = mySigma_i |
130 |
|
|
175 |
|
|
176 |
|
end subroutine createMixingList |
177 |
|
|
155 |
– |
|
178 |
|
subroutine do_lj_pair(atom1, atom2, d, rij, r2, pot, f, do_pot, do_stress) |
179 |
|
|
180 |
|
integer, intent(in) :: atom1, atom2 |
190 |
|
real( kind = dp ) :: pot_temp, dudr |
191 |
|
real( kind = dp ) :: sigma6 |
192 |
|
real( kind = dp ) :: epsilon |
171 |
– |
real( kind = dp ) :: rcut |
193 |
|
real( kind = dp ) :: r6 |
194 |
|
real( kind = dp ) :: t6 |
195 |
|
real( kind = dp ) :: t12 |
196 |
|
real( kind = dp ) :: delta |
197 |
|
|
198 |
< |
! Look up the correct parameters in the mixing matrix |
198 |
> |
|
199 |
> |
if (rij.lt.LJ_rcut) then |
200 |
> |
|
201 |
> |
! Look up the correct parameters in the mixing matrix |
202 |
|
#ifdef IS_MPI |
203 |
< |
sigma6 = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6 |
204 |
< |
epsilon = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon |
205 |
< |
delta = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta |
203 |
> |
sigma6 = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6 |
204 |
> |
epsilon = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon |
205 |
> |
delta = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta |
206 |
|
#else |
207 |
< |
sigma6 = ljMixed(atid(atom1),atid(atom2))%sigma6 |
208 |
< |
epsilon = ljMixed(atid(atom1),atid(atom2))%epsilon |
209 |
< |
delta = ljMixed(atid(atom1),atid(atom2))%delta |
207 |
> |
sigma6 = ljMixed(atid(atom1),atid(atom2))%sigma6 |
208 |
> |
epsilon = ljMixed(atid(atom1),atid(atom2))%epsilon |
209 |
> |
delta = ljMixed(atid(atom1),atid(atom2))%delta |
210 |
|
#endif |
187 |
– |
|
188 |
– |
call getRcut(rcut) |
189 |
– |
|
190 |
– |
r6 = r2 * r2 * r2 |
191 |
– |
|
192 |
– |
t6 = sigma6/ r6 |
193 |
– |
t12 = t6 * t6 |
194 |
– |
|
195 |
– |
if (rij.le.rcut) then |
211 |
|
|
212 |
+ |
r6 = r2 * r2 * r2 |
213 |
+ |
|
214 |
+ |
t6 = sigma6/ r6 |
215 |
+ |
t12 = t6 * t6 |
216 |
+ |
|
217 |
|
pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta |
218 |
< |
|
218 |
> |
|
219 |
|
dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij |
220 |
< |
|
220 |
> |
|
221 |
|
drdx = -d(1) / rij |
222 |
|
drdy = -d(2) / rij |
223 |
|
drdz = -d(3) / rij |
225 |
|
fx = dudr * drdx |
226 |
|
fy = dudr * drdy |
227 |
|
fz = dudr * drdz |
228 |
< |
|
228 |
> |
|
229 |
|
#ifdef IS_MPI |
230 |
|
if (do_pot) then |
231 |
|
pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5 |