45 |
|
|
46 |
|
!! @author Charles F. Vardeman II |
47 |
|
!! @author Matthew Meineke |
48 |
< |
!! @version $Id: doForces.F90,v 1.16 2005-05-17 02:09:06 kdaily Exp $, $Date: 2005-05-17 02:09:06 $, $Name: not supported by cvs2svn $, $Revision: 1.16 $ |
48 |
> |
!! @version $Id: doForces.F90,v 1.23 2005-07-03 20:53:43 chuckv Exp $, $Date: 2005-07-03 20:53:43 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $ |
49 |
|
|
50 |
|
|
51 |
|
module doForces |
73 |
|
|
74 |
|
#define __FORTRAN90 |
75 |
|
#include "UseTheForce/fSwitchingFunction.h" |
76 |
+ |
#include "UseTheForce/DarkSide/fInteractionMap.h" |
77 |
|
|
78 |
|
INTEGER, PARAMETER:: PREPAIR_LOOP = 1 |
79 |
|
INTEGER, PARAMETER:: PAIR_LOOP = 2 |
81 |
|
logical, save :: haveRlist = .false. |
82 |
|
logical, save :: haveNeighborList = .false. |
83 |
|
logical, save :: haveSIMvariables = .false. |
83 |
– |
logical, save :: havePropertyMap = .false. |
84 |
|
logical, save :: haveSaneForceField = .false. |
85 |
< |
|
85 |
> |
logical, save :: haveInteractionMap = .false. |
86 |
> |
|
87 |
|
logical, save :: FF_uses_DirectionalAtoms |
88 |
|
logical, save :: FF_uses_LennardJones |
89 |
|
logical, save :: FF_uses_Electrostatics |
90 |
|
logical, save :: FF_uses_Charges |
91 |
|
logical, save :: FF_uses_Dipoles |
92 |
|
logical, save :: FF_uses_Quadrupoles |
93 |
< |
logical, save :: FF_uses_sticky |
93 |
> |
logical, save :: FF_uses_Sticky |
94 |
> |
logical, save :: FF_uses_StickyPower |
95 |
|
logical, save :: FF_uses_GayBerne |
96 |
|
logical, save :: FF_uses_EAM |
97 |
|
logical, save :: FF_uses_Shapes |
105 |
|
logical, save :: SIM_uses_Dipoles |
106 |
|
logical, save :: SIM_uses_Quadrupoles |
107 |
|
logical, save :: SIM_uses_Sticky |
108 |
+ |
logical, save :: SIM_uses_StickyPower |
109 |
|
logical, save :: SIM_uses_GayBerne |
110 |
|
logical, save :: SIM_uses_EAM |
111 |
|
logical, save :: SIM_uses_Shapes |
116 |
|
logical, save :: SIM_uses_PBC |
117 |
|
logical, save :: SIM_uses_molecular_cutoffs |
118 |
|
|
119 |
< |
real(kind=dp), save :: rlist, rlistsq |
119 |
> |
!!!GO AWAY--------- |
120 |
> |
!!!!!real(kind=dp), save :: rlist, rlistsq |
121 |
|
|
122 |
|
public :: init_FF |
123 |
|
public :: do_force_loop |
124 |
< |
public :: setRlistDF |
124 |
> |
! public :: setRlistDF |
125 |
> |
!public :: addInteraction |
126 |
> |
!public :: setInteractionHash |
127 |
> |
!public :: getInteractionHash |
128 |
> |
public :: createInteractionMap |
129 |
> |
public :: createRcuts |
130 |
|
|
131 |
|
#ifdef PROFILE |
132 |
|
public :: getforcetime |
135 |
|
integer :: nLoops |
136 |
|
#endif |
137 |
|
|
138 |
< |
type :: Properties |
139 |
< |
logical :: is_Directional = .false. |
140 |
< |
logical :: is_LennardJones = .false. |
141 |
< |
logical :: is_Electrostatic = .false. |
142 |
< |
logical :: is_Charge = .false. |
143 |
< |
logical :: is_Dipole = .false. |
144 |
< |
logical :: is_Quadrupole = .false. |
145 |
< |
logical :: is_Sticky = .false. |
137 |
< |
logical :: is_GayBerne = .false. |
138 |
< |
logical :: is_EAM = .false. |
139 |
< |
logical :: is_Shape = .false. |
140 |
< |
logical :: is_FLARB = .false. |
141 |
< |
end type Properties |
138 |
> |
type, public :: Interaction |
139 |
> |
integer :: InteractionHash |
140 |
> |
real(kind=dp) :: rList = 0.0_dp |
141 |
> |
real(kind=dp) :: rListSq = 0.0_dp |
142 |
> |
end type Interaction |
143 |
> |
|
144 |
> |
type(Interaction), dimension(:,:),allocatable :: InteractionMap |
145 |
> |
|
146 |
|
|
147 |
< |
type(Properties), dimension(:),allocatable :: PropertyMap |
144 |
< |
|
147 |
> |
|
148 |
|
contains |
149 |
|
|
150 |
< |
subroutine setRlistDF( this_rlist ) |
151 |
< |
|
149 |
< |
real(kind=dp) :: this_rlist |
150 |
< |
|
151 |
< |
rlist = this_rlist |
152 |
< |
rlistsq = rlist * rlist |
153 |
< |
|
154 |
< |
haveRlist = .true. |
155 |
< |
|
156 |
< |
end subroutine setRlistDF |
157 |
< |
|
158 |
< |
subroutine createPropertyMap(status) |
150 |
> |
|
151 |
> |
subroutine createInteractionMap(status) |
152 |
|
integer :: nAtypes |
153 |
|
integer :: status |
154 |
|
integer :: i |
155 |
< |
logical :: thisProperty |
156 |
< |
real (kind=DP) :: thisDPproperty |
157 |
< |
|
158 |
< |
status = 0 |
159 |
< |
|
155 |
> |
integer :: j |
156 |
> |
integer :: ihash |
157 |
> |
real(kind=dp) :: myRcut |
158 |
> |
! Test Types |
159 |
> |
logical :: i_is_LJ |
160 |
> |
logical :: i_is_Elect |
161 |
> |
logical :: i_is_Sticky |
162 |
> |
logical :: i_is_StickyP |
163 |
> |
logical :: i_is_GB |
164 |
> |
logical :: i_is_EAM |
165 |
> |
logical :: i_is_Shape |
166 |
> |
logical :: j_is_LJ |
167 |
> |
logical :: j_is_Elect |
168 |
> |
logical :: j_is_Sticky |
169 |
> |
logical :: j_is_StickyP |
170 |
> |
logical :: j_is_GB |
171 |
> |
logical :: j_is_EAM |
172 |
> |
logical :: j_is_Shape |
173 |
> |
|
174 |
> |
|
175 |
> |
if (.not. associated(atypes)) then |
176 |
> |
call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!") |
177 |
> |
status = -1 |
178 |
> |
return |
179 |
> |
endif |
180 |
> |
|
181 |
|
nAtypes = getSize(atypes) |
182 |
< |
|
182 |
> |
|
183 |
|
if (nAtypes == 0) then |
184 |
|
status = -1 |
185 |
|
return |
186 |
|
end if |
173 |
– |
|
174 |
– |
if (.not. allocated(PropertyMap)) then |
175 |
– |
allocate(PropertyMap(nAtypes)) |
176 |
– |
endif |
187 |
|
|
188 |
+ |
if (.not. allocated(InteractionMap)) then |
189 |
+ |
allocate(InteractionMap(nAtypes,nAtypes)) |
190 |
+ |
endif |
191 |
+ |
|
192 |
|
do i = 1, nAtypes |
193 |
< |
call getElementProperty(atypes, i, "is_Directional", thisProperty) |
194 |
< |
PropertyMap(i)%is_Directional = thisProperty |
193 |
> |
call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) |
194 |
> |
call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect) |
195 |
> |
call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky) |
196 |
> |
call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP) |
197 |
> |
call getElementProperty(atypes, i, "is_GayBerne", i_is_GB) |
198 |
> |
call getElementProperty(atypes, i, "is_EAM", i_is_EAM) |
199 |
> |
call getElementProperty(atypes, i, "is_Shape", i_is_Shape) |
200 |
|
|
201 |
< |
call getElementProperty(atypes, i, "is_LennardJones", thisProperty) |
183 |
< |
PropertyMap(i)%is_LennardJones = thisProperty |
184 |
< |
|
185 |
< |
call getElementProperty(atypes, i, "is_Electrostatic", thisProperty) |
186 |
< |
PropertyMap(i)%is_Electrostatic = thisProperty |
201 |
> |
do j = i, nAtypes |
202 |
|
|
203 |
< |
call getElementProperty(atypes, i, "is_Charge", thisProperty) |
204 |
< |
PropertyMap(i)%is_Charge = thisProperty |
190 |
< |
|
191 |
< |
call getElementProperty(atypes, i, "is_Dipole", thisProperty) |
192 |
< |
PropertyMap(i)%is_Dipole = thisProperty |
203 |
> |
iHash = 0 |
204 |
> |
myRcut = 0.0_dp |
205 |
|
|
206 |
< |
call getElementProperty(atypes, i, "is_Quadrupole", thisProperty) |
207 |
< |
PropertyMap(i)%is_Quadrupole = thisProperty |
206 |
> |
call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ) |
207 |
> |
call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect) |
208 |
> |
call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky) |
209 |
> |
call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP) |
210 |
> |
call getElementProperty(atypes, j, "is_GayBerne", j_is_GB) |
211 |
> |
call getElementProperty(atypes, j, "is_EAM", j_is_EAM) |
212 |
> |
call getElementProperty(atypes, j, "is_Shape", j_is_Shape) |
213 |
|
|
214 |
< |
call getElementProperty(atypes, i, "is_Sticky", thisProperty) |
215 |
< |
PropertyMap(i)%is_Sticky = thisProperty |
214 |
> |
if (i_is_LJ .and. j_is_LJ) then |
215 |
> |
iHash = ior(iHash, LJ_PAIR) |
216 |
> |
endif |
217 |
> |
|
218 |
> |
if (i_is_Elect .and. j_is_Elect) then |
219 |
> |
iHash = ior(iHash, ELECTROSTATIC_PAIR) |
220 |
> |
endif |
221 |
> |
|
222 |
> |
if (i_is_Sticky .and. j_is_Sticky) then |
223 |
> |
iHash = ior(iHash, STICKY_PAIR) |
224 |
> |
endif |
225 |
|
|
226 |
< |
call getElementProperty(atypes, i, "is_GayBerne", thisProperty) |
227 |
< |
PropertyMap(i)%is_GayBerne = thisProperty |
226 |
> |
if (i_is_StickyP .and. j_is_StickyP) then |
227 |
> |
iHash = ior(iHash, STICKYPOWER_PAIR) |
228 |
> |
endif |
229 |
|
|
230 |
< |
call getElementProperty(atypes, i, "is_EAM", thisProperty) |
231 |
< |
PropertyMap(i)%is_EAM = thisProperty |
230 |
> |
if (i_is_EAM .and. j_is_EAM) then |
231 |
> |
iHash = ior(iHash, EAM_PAIR) |
232 |
> |
endif |
233 |
|
|
234 |
< |
call getElementProperty(atypes, i, "is_Shape", thisProperty) |
235 |
< |
PropertyMap(i)%is_Shape = thisProperty |
234 |
> |
if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR) |
235 |
> |
if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ) |
236 |
> |
if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ) |
237 |
|
|
238 |
< |
call getElementProperty(atypes, i, "is_FLARB", thisProperty) |
239 |
< |
PropertyMap(i)%is_FLARB = thisProperty |
238 |
> |
if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR) |
239 |
> |
if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ) |
240 |
> |
if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ) |
241 |
> |
|
242 |
> |
|
243 |
> |
InteractionMap(i,j)%InteractionHash = iHash |
244 |
> |
InteractionMap(j,i)%InteractionHash = iHash |
245 |
> |
|
246 |
> |
end do |
247 |
> |
|
248 |
|
end do |
249 |
+ |
end subroutine createInteractionMap |
250 |
|
|
251 |
< |
havePropertyMap = .true. |
251 |
> |
! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff. |
252 |
> |
subroutine createRcuts(defaultRList) |
253 |
> |
real(kind=dp), intent(in), optional :: defaultRList |
254 |
> |
integer :: iMap |
255 |
> |
integer :: map_i,map_j |
256 |
> |
real(kind=dp) :: thisRCut = 0.0_dp |
257 |
> |
real(kind=dp) :: actualCutoff = 0.0_dp |
258 |
> |
integer :: nAtypes |
259 |
|
|
260 |
< |
end subroutine createPropertyMap |
260 |
> |
if(.not. allocated(InteractionMap)) return |
261 |
|
|
262 |
+ |
nAtypes = getSize(atypes) |
263 |
+ |
! If we pass a default rcut, set all atypes to that cutoff distance |
264 |
+ |
if(present(defaultRList)) then |
265 |
+ |
InteractionMap(:,:)%rList = defaultRList |
266 |
+ |
InteractionMap(:,:)%rListSq = defaultRList*defaultRList |
267 |
+ |
haveRlist = .true. |
268 |
+ |
return |
269 |
+ |
end if |
270 |
+ |
|
271 |
+ |
do map_i = 1,nAtypes |
272 |
+ |
do map_j = map_i,nAtypes |
273 |
+ |
iMap = InteractionMap(map_i, map_j)%InteractionHash |
274 |
+ |
|
275 |
+ |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
276 |
+ |
! thisRCut = getLJCutOff(map_i,map_j) |
277 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
278 |
+ |
endif |
279 |
+ |
|
280 |
+ |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
281 |
+ |
! thisRCut = getElectrostaticCutOff(map_i,map_j) |
282 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
283 |
+ |
endif |
284 |
+ |
|
285 |
+ |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
286 |
+ |
! thisRCut = getStickyCutOff(map_i,map_j) |
287 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
288 |
+ |
endif |
289 |
+ |
|
290 |
+ |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
291 |
+ |
! thisRCut = getStickyPowerCutOff(map_i,map_j) |
292 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
293 |
+ |
endif |
294 |
+ |
|
295 |
+ |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
296 |
+ |
! thisRCut = getGayberneCutOff(map_i,map_j) |
297 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
298 |
+ |
endif |
299 |
+ |
|
300 |
+ |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
301 |
+ |
! thisRCut = getGaybrneLJCutOff(map_i,map_j) |
302 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
303 |
+ |
endif |
304 |
+ |
|
305 |
+ |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
306 |
+ |
! thisRCut = getEAMCutOff(map_i,map_j) |
307 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
308 |
+ |
endif |
309 |
+ |
|
310 |
+ |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
311 |
+ |
! thisRCut = getShapeCutOff(map_i,map_j) |
312 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
313 |
+ |
endif |
314 |
+ |
|
315 |
+ |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
316 |
+ |
! thisRCut = getShapeLJCutOff(map_i,map_j) |
317 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
318 |
+ |
endif |
319 |
+ |
InteractionMap(map_i, map_j)%rList = actualCutoff |
320 |
+ |
InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff |
321 |
+ |
end do |
322 |
+ |
end do |
323 |
+ |
haveRlist = .true. |
324 |
+ |
end subroutine createRcuts |
325 |
+ |
|
326 |
+ |
|
327 |
+ |
!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF |
328 |
+ |
!!$ subroutine setRlistDF( this_rlist ) |
329 |
+ |
!!$ |
330 |
+ |
!!$ real(kind=dp) :: this_rlist |
331 |
+ |
!!$ |
332 |
+ |
!!$ rlist = this_rlist |
333 |
+ |
!!$ rlistsq = rlist * rlist |
334 |
+ |
!!$ |
335 |
+ |
!!$ haveRlist = .true. |
336 |
+ |
!!$ |
337 |
+ |
!!$ end subroutine setRlistDF |
338 |
+ |
|
339 |
+ |
|
340 |
|
subroutine setSimVariables() |
341 |
|
SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms() |
342 |
|
SIM_uses_LennardJones = SimUsesLennardJones() |
344 |
|
SIM_uses_Charges = SimUsesCharges() |
345 |
|
SIM_uses_Dipoles = SimUsesDipoles() |
346 |
|
SIM_uses_Sticky = SimUsesSticky() |
347 |
+ |
SIM_uses_StickyPower = SimUsesStickyPower() |
348 |
|
SIM_uses_GayBerne = SimUsesGayBerne() |
349 |
|
SIM_uses_EAM = SimUsesEAM() |
350 |
|
SIM_uses_Shapes = SimUsesShapes() |
365 |
|
integer :: myStatus |
366 |
|
|
367 |
|
error = 0 |
244 |
– |
|
245 |
– |
if (.not. havePropertyMap) then |
368 |
|
|
369 |
+ |
if (.not. haveInteractionMap) then |
370 |
+ |
|
371 |
|
myStatus = 0 |
372 |
|
|
373 |
< |
call createPropertyMap(myStatus) |
373 |
> |
call createInteractionMap(myStatus) |
374 |
|
|
375 |
|
if (myStatus .ne. 0) then |
376 |
< |
write(default_error, *) 'createPropertyMap failed in doForces!' |
376 |
> |
write(default_error, *) 'createInteractionMap failed in doForces!' |
377 |
|
error = -1 |
378 |
|
return |
379 |
|
endif |
410 |
|
#endif |
411 |
|
return |
412 |
|
end subroutine doReadyCheck |
289 |
– |
|
413 |
|
|
414 |
+ |
|
415 |
|
subroutine init_FF(use_RF_c, thisStat) |
416 |
|
|
417 |
|
logical, intent(in) :: use_RF_c |
426 |
|
|
427 |
|
!! Fortran's version of a cast: |
428 |
|
FF_uses_RF = use_RF_c |
429 |
< |
|
429 |
> |
|
430 |
|
!! init_FF is called *after* all of the atom types have been |
431 |
|
!! defined in atype_module using the new_atype subroutine. |
432 |
|
!! |
433 |
|
!! this will scan through the known atypes and figure out what |
434 |
|
!! interactions are used by the force field. |
435 |
< |
|
435 |
> |
|
436 |
|
FF_uses_DirectionalAtoms = .false. |
437 |
|
FF_uses_LennardJones = .false. |
438 |
|
FF_uses_Electrostatics = .false. |
439 |
|
FF_uses_Charges = .false. |
440 |
|
FF_uses_Dipoles = .false. |
441 |
|
FF_uses_Sticky = .false. |
442 |
+ |
FF_uses_StickyPower = .false. |
443 |
|
FF_uses_GayBerne = .false. |
444 |
|
FF_uses_EAM = .false. |
445 |
|
FF_uses_Shapes = .false. |
446 |
|
FF_uses_FLARB = .false. |
447 |
< |
|
447 |
> |
|
448 |
|
call getMatchingElementList(atypes, "is_Directional", .true., & |
449 |
|
nMatches, MatchList) |
450 |
|
if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true. |
452 |
|
call getMatchingElementList(atypes, "is_LennardJones", .true., & |
453 |
|
nMatches, MatchList) |
454 |
|
if (nMatches .gt. 0) FF_uses_LennardJones = .true. |
455 |
< |
|
455 |
> |
|
456 |
|
call getMatchingElementList(atypes, "is_Electrostatic", .true., & |
457 |
|
nMatches, MatchList) |
458 |
|
if (nMatches .gt. 0) then |
465 |
|
FF_uses_Charges = .true. |
466 |
|
FF_uses_Electrostatics = .true. |
467 |
|
endif |
468 |
< |
|
468 |
> |
|
469 |
|
call getMatchingElementList(atypes, "is_Dipole", .true., & |
470 |
|
nMatches, MatchList) |
471 |
|
if (nMatches .gt. 0) then |
481 |
|
FF_uses_Electrostatics = .true. |
482 |
|
FF_uses_DirectionalAtoms = .true. |
483 |
|
endif |
484 |
< |
|
484 |
> |
|
485 |
|
call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, & |
486 |
|
MatchList) |
487 |
|
if (nMatches .gt. 0) then |
488 |
|
FF_uses_Sticky = .true. |
489 |
|
FF_uses_DirectionalAtoms = .true. |
490 |
|
endif |
491 |
+ |
|
492 |
+ |
call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, & |
493 |
+ |
MatchList) |
494 |
+ |
if (nMatches .gt. 0) then |
495 |
+ |
FF_uses_StickyPower = .true. |
496 |
+ |
FF_uses_DirectionalAtoms = .true. |
497 |
+ |
endif |
498 |
|
|
499 |
|
call getMatchingElementList(atypes, "is_GayBerne", .true., & |
500 |
|
nMatches, MatchList) |
502 |
|
FF_uses_GayBerne = .true. |
503 |
|
FF_uses_DirectionalAtoms = .true. |
504 |
|
endif |
505 |
< |
|
505 |
> |
|
506 |
|
call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList) |
507 |
|
if (nMatches .gt. 0) FF_uses_EAM = .true. |
508 |
< |
|
508 |
> |
|
509 |
|
call getMatchingElementList(atypes, "is_Shape", .true., & |
510 |
|
nMatches, MatchList) |
511 |
|
if (nMatches .gt. 0) then |
519 |
|
|
520 |
|
!! Assume sanity (for the sake of argument) |
521 |
|
haveSaneForceField = .true. |
522 |
< |
|
522 |
> |
|
523 |
|
!! check to make sure the FF_uses_RF setting makes sense |
524 |
< |
|
524 |
> |
|
525 |
|
if (FF_uses_dipoles) then |
526 |
|
if (FF_uses_RF) then |
527 |
|
dielect = getDielect() |
534 |
|
haveSaneForceField = .false. |
535 |
|
return |
536 |
|
endif |
537 |
< |
endif |
537 |
> |
endif |
538 |
|
|
539 |
|
!sticky module does not contain check_sticky_FF anymore |
540 |
|
!if (FF_uses_sticky) then |
547 |
|
!endif |
548 |
|
|
549 |
|
if (FF_uses_EAM) then |
550 |
< |
call init_EAM_FF(my_status) |
550 |
> |
call init_EAM_FF(my_status) |
551 |
|
if (my_status /= 0) then |
552 |
|
write(default_error, *) "init_EAM_FF returned a bad status" |
553 |
|
thisStat = -1 |
567 |
|
|
568 |
|
if (FF_uses_GayBerne .and. FF_uses_LennardJones) then |
569 |
|
endif |
570 |
< |
|
570 |
> |
|
571 |
|
if (.not. haveNeighborList) then |
572 |
|
!! Create neighbor lists |
573 |
|
call expandNeighborList(nLocal, my_status) |
577 |
|
return |
578 |
|
endif |
579 |
|
haveNeighborList = .true. |
580 |
< |
endif |
581 |
< |
|
580 |
> |
endif |
581 |
> |
|
582 |
|
end subroutine init_FF |
451 |
– |
|
583 |
|
|
584 |
+ |
|
585 |
|
!! Does force loop over i,j pairs. Calls do_pair to calculates forces. |
586 |
|
!-------------------------------------------------------------> |
587 |
|
subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, & |
631 |
|
integer :: localError |
632 |
|
integer :: propPack_i, propPack_j |
633 |
|
integer :: loopStart, loopEnd, loop |
634 |
< |
|
634 |
> |
integer :: iMap |
635 |
|
real(kind=dp) :: listSkin = 1.0 |
636 |
< |
|
636 |
> |
|
637 |
|
!! initialize local variables |
638 |
< |
|
638 |
> |
|
639 |
|
#ifdef IS_MPI |
640 |
|
pot_local = 0.0_dp |
641 |
|
nAtomsInRow = getNatomsInRow(plan_atom_row) |
645 |
|
#else |
646 |
|
natoms = nlocal |
647 |
|
#endif |
648 |
< |
|
648 |
> |
|
649 |
|
call doReadyCheck(localError) |
650 |
|
if ( localError .ne. 0 ) then |
651 |
|
call handleError("do_force_loop", "Not Initialized") |
653 |
|
return |
654 |
|
end if |
655 |
|
call zero_work_arrays() |
656 |
< |
|
656 |
> |
|
657 |
|
do_pot = do_pot_c |
658 |
|
do_stress = do_stress_c |
659 |
< |
|
659 |
> |
|
660 |
|
! Gather all information needed by all force loops: |
661 |
< |
|
661 |
> |
|
662 |
|
#ifdef IS_MPI |
663 |
< |
|
663 |
> |
|
664 |
|
call gather(q, q_Row, plan_atom_row_3d) |
665 |
|
call gather(q, q_Col, plan_atom_col_3d) |
666 |
|
|
667 |
|
call gather(q_group, q_group_Row, plan_group_row_3d) |
668 |
|
call gather(q_group, q_group_Col, plan_group_col_3d) |
669 |
< |
|
669 |
> |
|
670 |
|
if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then |
671 |
|
call gather(eFrame, eFrame_Row, plan_atom_row_rotation) |
672 |
|
call gather(eFrame, eFrame_Col, plan_atom_col_rotation) |
673 |
< |
|
673 |
> |
|
674 |
|
call gather(A, A_Row, plan_atom_row_rotation) |
675 |
|
call gather(A, A_Col, plan_atom_col_rotation) |
676 |
|
endif |
677 |
< |
|
677 |
> |
|
678 |
|
#endif |
679 |
< |
|
679 |
> |
|
680 |
|
!! Begin force loop timing: |
681 |
|
#ifdef PROFILE |
682 |
|
call cpu_time(forceTimeInitial) |
683 |
|
nloops = nloops + 1 |
684 |
|
#endif |
685 |
< |
|
685 |
> |
|
686 |
|
loopEnd = PAIR_LOOP |
687 |
|
if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then |
688 |
|
loopStart = PREPAIR_LOOP |
697 |
|
if (loop .eq. loopStart) then |
698 |
|
#ifdef IS_MPI |
699 |
|
call checkNeighborList(nGroupsInRow, q_group_row, listSkin, & |
700 |
< |
update_nlist) |
700 |
> |
update_nlist) |
701 |
|
#else |
702 |
|
call checkNeighborList(nGroups, q_group, listSkin, & |
703 |
< |
update_nlist) |
703 |
> |
update_nlist) |
704 |
|
#endif |
705 |
|
endif |
706 |
< |
|
706 |
> |
|
707 |
|
if (update_nlist) then |
708 |
|
!! save current configuration and construct neighbor list |
709 |
|
#ifdef IS_MPI |
714 |
|
neighborListSize = size(list) |
715 |
|
nlist = 0 |
716 |
|
endif |
717 |
< |
|
717 |
> |
|
718 |
|
istart = 1 |
719 |
|
#ifdef IS_MPI |
720 |
|
iend = nGroupsInRow |
724 |
|
outer: do i = istart, iend |
725 |
|
|
726 |
|
if (update_nlist) point(i) = nlist + 1 |
727 |
< |
|
727 |
> |
|
728 |
|
n_in_i = groupStartRow(i+1) - groupStartRow(i) |
729 |
< |
|
729 |
> |
|
730 |
|
if (update_nlist) then |
731 |
|
#ifdef IS_MPI |
732 |
|
jstart = 1 |
741 |
|
! make sure group i has neighbors |
742 |
|
if (jstart .gt. jend) cycle outer |
743 |
|
endif |
744 |
< |
|
744 |
> |
|
745 |
|
do jnab = jstart, jend |
746 |
|
if (update_nlist) then |
747 |
|
j = jnab |
757 |
|
q_group(:,j), d_grp, rgrpsq) |
758 |
|
#endif |
759 |
|
|
760 |
< |
if (rgrpsq < rlistsq) then |
760 |
> |
if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then |
761 |
|
if (update_nlist) then |
762 |
|
nlist = nlist + 1 |
763 |
< |
|
763 |
> |
|
764 |
|
if (nlist > neighborListSize) then |
765 |
|
#ifdef IS_MPI |
766 |
|
call expandNeighborList(nGroupsInRow, listerror) |
774 |
|
end if |
775 |
|
neighborListSize = size(list) |
776 |
|
endif |
777 |
< |
|
777 |
> |
|
778 |
|
list(nlist) = j |
779 |
|
endif |
780 |
< |
|
780 |
> |
|
781 |
|
if (loop .eq. PAIR_LOOP) then |
782 |
|
vij = 0.0d0 |
783 |
|
fij(1:3) = 0.0d0 |
784 |
|
endif |
785 |
< |
|
785 |
> |
|
786 |
|
call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, & |
787 |
|
in_switching_region) |
788 |
< |
|
788 |
> |
|
789 |
|
n_in_j = groupStartCol(j+1) - groupStartCol(j) |
790 |
< |
|
790 |
> |
|
791 |
|
do ia = groupStartRow(i), groupStartRow(i+1)-1 |
792 |
< |
|
792 |
> |
|
793 |
|
atom1 = groupListRow(ia) |
794 |
< |
|
794 |
> |
|
795 |
|
inner: do jb = groupStartCol(j), groupStartCol(j+1)-1 |
796 |
< |
|
796 |
> |
|
797 |
|
atom2 = groupListCol(jb) |
798 |
< |
|
798 |
> |
|
799 |
|
if (skipThisPair(atom1, atom2)) cycle inner |
800 |
|
|
801 |
|
if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then |
837 |
|
endif |
838 |
|
enddo inner |
839 |
|
enddo |
840 |
< |
|
840 |
> |
|
841 |
|
if (loop .eq. PAIR_LOOP) then |
842 |
|
if (in_switching_region) then |
843 |
|
swderiv = vij*dswdr/rgrp |
844 |
|
fij(1) = fij(1) + swderiv*d_grp(1) |
845 |
|
fij(2) = fij(2) + swderiv*d_grp(2) |
846 |
|
fij(3) = fij(3) + swderiv*d_grp(3) |
847 |
< |
|
847 |
> |
|
848 |
|
do ia=groupStartRow(i), groupStartRow(i+1)-1 |
849 |
|
atom1=groupListRow(ia) |
850 |
|
mf = mfactRow(atom1) |
858 |
|
f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf |
859 |
|
#endif |
860 |
|
enddo |
861 |
< |
|
861 |
> |
|
862 |
|
do jb=groupStartCol(j), groupStartCol(j+1)-1 |
863 |
|
atom2=groupListCol(jb) |
864 |
|
mf = mfactCol(atom2) |
873 |
|
#endif |
874 |
|
enddo |
875 |
|
endif |
876 |
< |
|
876 |
> |
|
877 |
|
if (do_stress) call add_stress_tensor(d_grp, fij) |
878 |
|
endif |
879 |
|
end if |
880 |
|
enddo |
881 |
|
enddo outer |
882 |
< |
|
882 |
> |
|
883 |
|
if (update_nlist) then |
884 |
|
#ifdef IS_MPI |
885 |
|
point(nGroupsInRow + 1) = nlist + 1 |
893 |
|
update_nlist = .false. |
894 |
|
endif |
895 |
|
endif |
896 |
< |
|
896 |
> |
|
897 |
|
if (loop .eq. PREPAIR_LOOP) then |
898 |
|
call do_preforce(nlocal, pot) |
899 |
|
endif |
900 |
< |
|
900 |
> |
|
901 |
|
enddo |
902 |
< |
|
902 |
> |
|
903 |
|
!! Do timing |
904 |
|
#ifdef PROFILE |
905 |
|
call cpu_time(forceTimeFinal) |
906 |
|
forceTime = forceTime + forceTimeFinal - forceTimeInitial |
907 |
|
#endif |
908 |
< |
|
908 |
> |
|
909 |
|
#ifdef IS_MPI |
910 |
|
!!distribute forces |
911 |
< |
|
911 |
> |
|
912 |
|
f_temp = 0.0_dp |
913 |
|
call scatter(f_Row,f_temp,plan_atom_row_3d) |
914 |
|
do i = 1,nlocal |
915 |
|
f(1:3,i) = f(1:3,i) + f_temp(1:3,i) |
916 |
|
end do |
917 |
< |
|
917 |
> |
|
918 |
|
f_temp = 0.0_dp |
919 |
|
call scatter(f_Col,f_temp,plan_atom_col_3d) |
920 |
|
do i = 1,nlocal |
921 |
|
f(1:3,i) = f(1:3,i) + f_temp(1:3,i) |
922 |
|
end do |
923 |
< |
|
923 |
> |
|
924 |
|
if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then |
925 |
|
t_temp = 0.0_dp |
926 |
|
call scatter(t_Row,t_temp,plan_atom_row_3d) |
929 |
|
end do |
930 |
|
t_temp = 0.0_dp |
931 |
|
call scatter(t_Col,t_temp,plan_atom_col_3d) |
932 |
< |
|
932 |
> |
|
933 |
|
do i = 1,nlocal |
934 |
|
t(1:3,i) = t(1:3,i) + t_temp(1:3,i) |
935 |
|
end do |
936 |
|
endif |
937 |
< |
|
937 |
> |
|
938 |
|
if (do_pot) then |
939 |
|
! scatter/gather pot_row into the members of my column |
940 |
|
call scatter(pot_Row, pot_Temp, plan_atom_row) |
941 |
< |
|
941 |
> |
|
942 |
|
! scatter/gather pot_local into all other procs |
943 |
|
! add resultant to get total pot |
944 |
|
do i = 1, nlocal |
945 |
|
pot_local = pot_local + pot_Temp(i) |
946 |
|
enddo |
947 |
< |
|
947 |
> |
|
948 |
|
pot_Temp = 0.0_DP |
949 |
< |
|
949 |
> |
|
950 |
|
call scatter(pot_Col, pot_Temp, plan_atom_col) |
951 |
|
do i = 1, nlocal |
952 |
|
pot_local = pot_local + pot_Temp(i) |
953 |
|
enddo |
954 |
< |
|
954 |
> |
|
955 |
|
endif |
956 |
|
#endif |
957 |
< |
|
957 |
> |
|
958 |
|
if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then |
959 |
< |
|
959 |
> |
|
960 |
|
if (FF_uses_RF .and. SIM_uses_RF) then |
961 |
< |
|
961 |
> |
|
962 |
|
#ifdef IS_MPI |
963 |
|
call scatter(rf_Row,rf,plan_atom_row_3d) |
964 |
|
call scatter(rf_Col,rf_Temp,plan_atom_col_3d) |
966 |
|
rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i) |
967 |
|
end do |
968 |
|
#endif |
969 |
< |
|
969 |
> |
|
970 |
|
do i = 1, nLocal |
971 |
< |
|
971 |
> |
|
972 |
|
rfpot = 0.0_DP |
973 |
|
#ifdef IS_MPI |
974 |
|
me_i = atid_row(i) |
975 |
|
#else |
976 |
|
me_i = atid(i) |
977 |
|
#endif |
978 |
+ |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
979 |
|
|
980 |
< |
if (PropertyMap(me_i)%is_Dipole) then |
981 |
< |
|
980 |
> |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
981 |
> |
|
982 |
|
mu_i = getDipoleMoment(me_i) |
983 |
< |
|
983 |
> |
|
984 |
|
!! The reaction field needs to include a self contribution |
985 |
|
!! to the field: |
986 |
|
call accumulate_self_rf(i, mu_i, eFrame) |
991 |
|
pot_local = pot_local + rfpot |
992 |
|
#else |
993 |
|
pot = pot + rfpot |
994 |
< |
|
994 |
> |
|
995 |
|
#endif |
996 |
< |
endif |
996 |
> |
endif |
997 |
|
enddo |
998 |
|
endif |
999 |
|
endif |
1000 |
< |
|
1001 |
< |
|
1000 |
> |
|
1001 |
> |
|
1002 |
|
#ifdef IS_MPI |
1003 |
< |
|
1003 |
> |
|
1004 |
|
if (do_pot) then |
1005 |
|
pot = pot + pot_local |
1006 |
|
!! we assume the c code will do the allreduce to get the total potential |
1007 |
|
!! we could do it right here if we needed to... |
1008 |
|
endif |
1009 |
< |
|
1009 |
> |
|
1010 |
|
if (do_stress) then |
1011 |
|
call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, & |
1012 |
|
mpi_comm_world,mpi_err) |
1013 |
|
call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, & |
1014 |
|
mpi_comm_world,mpi_err) |
1015 |
|
endif |
1016 |
< |
|
1016 |
> |
|
1017 |
|
#else |
1018 |
< |
|
1018 |
> |
|
1019 |
|
if (do_stress) then |
1020 |
|
tau = tau_Temp |
1021 |
|
virial = virial_Temp |
1022 |
|
endif |
1023 |
< |
|
1023 |
> |
|
1024 |
|
#endif |
1025 |
< |
|
1025 |
> |
|
1026 |
|
end subroutine do_force_loop |
1027 |
< |
|
1027 |
> |
|
1028 |
|
subroutine do_pair(i, j, rijsq, d, sw, do_pot, & |
1029 |
|
eFrame, A, f, t, pot, vpair, fpair) |
1030 |
|
|
1041 |
|
real ( kind = dp ), intent(inout) :: rijsq |
1042 |
|
real ( kind = dp ) :: r |
1043 |
|
real ( kind = dp ), intent(inout) :: d(3) |
1044 |
+ |
real ( kind = dp ) :: ebalance |
1045 |
|
integer :: me_i, me_j |
1046 |
|
|
1047 |
+ |
integer :: iMap |
1048 |
+ |
|
1049 |
|
r = sqrt(rijsq) |
1050 |
|
vpair = 0.0d0 |
1051 |
|
fpair(1:3) = 0.0d0 |
1058 |
|
me_j = atid(j) |
1059 |
|
#endif |
1060 |
|
|
1061 |
< |
! write(*,*) i, j, me_i, me_j |
1062 |
< |
|
1063 |
< |
if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then |
1064 |
< |
|
929 |
< |
if ( PropertyMap(me_i)%is_LennardJones .and. & |
930 |
< |
PropertyMap(me_j)%is_LennardJones ) then |
931 |
< |
call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) |
932 |
< |
endif |
933 |
< |
|
1061 |
> |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
1062 |
> |
|
1063 |
> |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
1064 |
> |
call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) |
1065 |
|
endif |
935 |
– |
|
936 |
– |
if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then |
937 |
– |
|
938 |
– |
if (PropertyMap(me_i)%is_Electrostatic .and. & |
939 |
– |
PropertyMap(me_j)%is_Electrostatic) then |
940 |
– |
call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & |
941 |
– |
pot, eFrame, f, t, do_pot) |
942 |
– |
endif |
943 |
– |
|
944 |
– |
if (FF_uses_dipoles .and. SIM_uses_dipoles) then |
945 |
– |
if ( PropertyMap(me_i)%is_Dipole .and. & |
946 |
– |
PropertyMap(me_j)%is_Dipole) then |
947 |
– |
if (FF_uses_RF .and. SIM_uses_RF) then |
948 |
– |
call accumulate_rf(i, j, r, eFrame, sw) |
949 |
– |
call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair) |
950 |
– |
endif |
951 |
– |
endif |
952 |
– |
endif |
953 |
– |
endif |
1066 |
|
|
1067 |
+ |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
1068 |
+ |
call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1069 |
+ |
pot, eFrame, f, t, do_pot) |
1070 |
|
|
1071 |
< |
if (FF_uses_Sticky .and. SIM_uses_sticky) then |
1071 |
> |
if (FF_uses_RF .and. SIM_uses_RF) then |
1072 |
|
|
1073 |
< |
if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then |
1074 |
< |
call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1075 |
< |
pot, A, f, t, do_pot) |
1073 |
> |
! CHECK ME (RF needs to know about all electrostatic types) |
1074 |
> |
call accumulate_rf(i, j, r, eFrame, sw) |
1075 |
> |
call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair) |
1076 |
|
endif |
1077 |
< |
|
1077 |
> |
|
1078 |
|
endif |
1079 |
|
|
1080 |
+ |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
1081 |
+ |
call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1082 |
+ |
pot, A, f, t, do_pot) |
1083 |
+ |
endif |
1084 |
|
|
1085 |
< |
if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then |
1086 |
< |
|
1087 |
< |
if ( PropertyMap(me_i)%is_GayBerne .and. & |
969 |
< |
PropertyMap(me_j)%is_GayBerne) then |
970 |
< |
call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
971 |
< |
pot, A, f, t, do_pot) |
972 |
< |
endif |
973 |
< |
|
1085 |
> |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
1086 |
> |
call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1087 |
> |
pot, A, f, t, do_pot) |
1088 |
|
endif |
1089 |
+ |
|
1090 |
+ |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
1091 |
+ |
call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1092 |
+ |
pot, A, f, t, do_pot) |
1093 |
+ |
endif |
1094 |
|
|
1095 |
< |
if (FF_uses_EAM .and. SIM_uses_EAM) then |
1096 |
< |
|
1097 |
< |
if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then |
979 |
< |
call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, & |
980 |
< |
do_pot) |
981 |
< |
endif |
982 |
< |
|
1095 |
> |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
1096 |
> |
! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1097 |
> |
! pot, A, f, t, do_pot) |
1098 |
|
endif |
1099 |
|
|
1100 |
+ |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
1101 |
+ |
call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, & |
1102 |
+ |
do_pot) |
1103 |
+ |
endif |
1104 |
|
|
1105 |
< |
! write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape |
1105 |
> |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
1106 |
> |
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1107 |
> |
pot, A, f, t, do_pot) |
1108 |
> |
endif |
1109 |
|
|
1110 |
< |
if (FF_uses_Shapes .and. SIM_uses_Shapes) then |
1111 |
< |
if ( PropertyMap(me_i)%is_Shape .and. & |
1112 |
< |
PropertyMap(me_j)%is_Shape ) then |
991 |
< |
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
992 |
< |
pot, A, f, t, do_pot) |
993 |
< |
endif |
994 |
< |
|
1110 |
> |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
1111 |
> |
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1112 |
> |
pot, A, f, t, do_pot) |
1113 |
|
endif |
1114 |
|
|
1115 |
|
end subroutine do_pair |
1117 |
|
subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, & |
1118 |
|
do_pot, do_stress, eFrame, A, f, t, pot) |
1119 |
|
|
1120 |
< |
real( kind = dp ) :: pot, sw |
1121 |
< |
real( kind = dp ), dimension(9,nLocal) :: eFrame |
1122 |
< |
real (kind=dp), dimension(9,nLocal) :: A |
1123 |
< |
real (kind=dp), dimension(3,nLocal) :: f |
1124 |
< |
real (kind=dp), dimension(3,nLocal) :: t |
1007 |
< |
|
1008 |
< |
logical, intent(inout) :: do_pot, do_stress |
1009 |
< |
integer, intent(in) :: i, j |
1010 |
< |
real ( kind = dp ), intent(inout) :: rijsq, rcijsq |
1011 |
< |
real ( kind = dp ) :: r, rc |
1012 |
< |
real ( kind = dp ), intent(inout) :: d(3), dc(3) |
1013 |
< |
|
1014 |
< |
logical :: is_EAM_i, is_EAM_j |
1015 |
< |
|
1016 |
< |
integer :: me_i, me_j |
1017 |
< |
|
1120 |
> |
real( kind = dp ) :: pot, sw |
1121 |
> |
real( kind = dp ), dimension(9,nLocal) :: eFrame |
1122 |
> |
real (kind=dp), dimension(9,nLocal) :: A |
1123 |
> |
real (kind=dp), dimension(3,nLocal) :: f |
1124 |
> |
real (kind=dp), dimension(3,nLocal) :: t |
1125 |
|
|
1126 |
< |
r = sqrt(rijsq) |
1127 |
< |
if (SIM_uses_molecular_cutoffs) then |
1128 |
< |
rc = sqrt(rcijsq) |
1129 |
< |
else |
1130 |
< |
rc = r |
1024 |
< |
endif |
1025 |
< |
|
1126 |
> |
logical, intent(inout) :: do_pot, do_stress |
1127 |
> |
integer, intent(in) :: i, j |
1128 |
> |
real ( kind = dp ), intent(inout) :: rijsq, rcijsq |
1129 |
> |
real ( kind = dp ) :: r, rc |
1130 |
> |
real ( kind = dp ), intent(inout) :: d(3), dc(3) |
1131 |
|
|
1132 |
+ |
integer :: me_i, me_j, iMap |
1133 |
+ |
|
1134 |
|
#ifdef IS_MPI |
1135 |
< |
me_i = atid_row(i) |
1136 |
< |
me_j = atid_col(j) |
1135 |
> |
me_i = atid_row(i) |
1136 |
> |
me_j = atid_col(j) |
1137 |
|
#else |
1138 |
< |
me_i = atid(i) |
1139 |
< |
me_j = atid(j) |
1138 |
> |
me_i = atid(i) |
1139 |
> |
me_j = atid(j) |
1140 |
|
#endif |
1034 |
– |
|
1035 |
– |
if (FF_uses_EAM .and. SIM_uses_EAM) then |
1036 |
– |
|
1037 |
– |
if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) & |
1038 |
– |
call calc_EAM_prepair_rho(i, j, d, r, rijsq ) |
1039 |
– |
|
1040 |
– |
endif |
1041 |
– |
|
1042 |
– |
end subroutine do_prepair |
1043 |
– |
|
1044 |
– |
|
1045 |
– |
subroutine do_preforce(nlocal,pot) |
1046 |
– |
integer :: nlocal |
1047 |
– |
real( kind = dp ) :: pot |
1048 |
– |
|
1049 |
– |
if (FF_uses_EAM .and. SIM_uses_EAM) then |
1050 |
– |
call calc_EAM_preforce_Frho(nlocal,pot) |
1051 |
– |
endif |
1052 |
– |
|
1053 |
– |
|
1054 |
– |
end subroutine do_preforce |
1055 |
– |
|
1056 |
– |
|
1057 |
– |
subroutine get_interatomic_vector(q_i, q_j, d, r_sq) |
1058 |
– |
|
1059 |
– |
real (kind = dp), dimension(3) :: q_i |
1060 |
– |
real (kind = dp), dimension(3) :: q_j |
1061 |
– |
real ( kind = dp ), intent(out) :: r_sq |
1062 |
– |
real( kind = dp ) :: d(3), scaled(3) |
1063 |
– |
integer i |
1064 |
– |
|
1065 |
– |
d(1:3) = q_j(1:3) - q_i(1:3) |
1066 |
– |
|
1067 |
– |
! Wrap back into periodic box if necessary |
1068 |
– |
if ( SIM_uses_PBC ) then |
1069 |
– |
|
1070 |
– |
if( .not.boxIsOrthorhombic ) then |
1071 |
– |
! calc the scaled coordinates. |
1072 |
– |
|
1073 |
– |
scaled = matmul(HmatInv, d) |
1074 |
– |
|
1075 |
– |
! wrap the scaled coordinates |
1076 |
– |
|
1077 |
– |
scaled = scaled - anint(scaled) |
1078 |
– |
|
1079 |
– |
|
1080 |
– |
! calc the wrapped real coordinates from the wrapped scaled |
1081 |
– |
! coordinates |
1082 |
– |
|
1083 |
– |
d = matmul(Hmat,scaled) |
1084 |
– |
|
1085 |
– |
else |
1086 |
– |
! calc the scaled coordinates. |
1087 |
– |
|
1088 |
– |
do i = 1, 3 |
1089 |
– |
scaled(i) = d(i) * HmatInv(i,i) |
1090 |
– |
|
1091 |
– |
! wrap the scaled coordinates |
1092 |
– |
|
1093 |
– |
scaled(i) = scaled(i) - anint(scaled(i)) |
1094 |
– |
|
1095 |
– |
! calc the wrapped real coordinates from the wrapped scaled |
1096 |
– |
! coordinates |
1097 |
– |
|
1098 |
– |
d(i) = scaled(i)*Hmat(i,i) |
1099 |
– |
enddo |
1100 |
– |
endif |
1101 |
– |
|
1102 |
– |
endif |
1103 |
– |
|
1104 |
– |
r_sq = dot_product(d,d) |
1105 |
– |
|
1106 |
– |
end subroutine get_interatomic_vector |
1107 |
– |
|
1108 |
– |
subroutine zero_work_arrays() |
1109 |
– |
|
1110 |
– |
#ifdef IS_MPI |
1111 |
– |
|
1112 |
– |
q_Row = 0.0_dp |
1113 |
– |
q_Col = 0.0_dp |
1141 |
|
|
1142 |
< |
q_group_Row = 0.0_dp |
1143 |
< |
q_group_Col = 0.0_dp |
1144 |
< |
|
1145 |
< |
eFrame_Row = 0.0_dp |
1146 |
< |
eFrame_Col = 0.0_dp |
1147 |
< |
|
1148 |
< |
A_Row = 0.0_dp |
1149 |
< |
A_Col = 0.0_dp |
1150 |
< |
|
1151 |
< |
f_Row = 0.0_dp |
1152 |
< |
f_Col = 0.0_dp |
1153 |
< |
f_Temp = 0.0_dp |
1154 |
< |
|
1155 |
< |
t_Row = 0.0_dp |
1156 |
< |
t_Col = 0.0_dp |
1157 |
< |
t_Temp = 0.0_dp |
1158 |
< |
|
1159 |
< |
pot_Row = 0.0_dp |
1160 |
< |
pot_Col = 0.0_dp |
1161 |
< |
pot_Temp = 0.0_dp |
1162 |
< |
|
1163 |
< |
rf_Row = 0.0_dp |
1164 |
< |
rf_Col = 0.0_dp |
1165 |
< |
rf_Temp = 0.0_dp |
1166 |
< |
|
1167 |
< |
#endif |
1168 |
< |
|
1169 |
< |
if (FF_uses_EAM .and. SIM_uses_EAM) then |
1170 |
< |
call clean_EAM() |
1171 |
< |
endif |
1172 |
< |
|
1173 |
< |
rf = 0.0_dp |
1174 |
< |
tau_Temp = 0.0_dp |
1175 |
< |
virial_Temp = 0.0_dp |
1176 |
< |
end subroutine zero_work_arrays |
1177 |
< |
|
1178 |
< |
function skipThisPair(atom1, atom2) result(skip_it) |
1179 |
< |
integer, intent(in) :: atom1 |
1180 |
< |
integer, intent(in), optional :: atom2 |
1181 |
< |
logical :: skip_it |
1182 |
< |
integer :: unique_id_1, unique_id_2 |
1183 |
< |
integer :: me_i,me_j |
1184 |
< |
integer :: i |
1185 |
< |
|
1186 |
< |
skip_it = .false. |
1187 |
< |
|
1188 |
< |
!! there are a number of reasons to skip a pair or a particle |
1189 |
< |
!! mostly we do this to exclude atoms who are involved in short |
1190 |
< |
!! range interactions (bonds, bends, torsions), but we also need |
1191 |
< |
!! to exclude some overcounted interactions that result from |
1192 |
< |
!! the parallel decomposition |
1193 |
< |
|
1142 |
> |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
1143 |
> |
|
1144 |
> |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
1145 |
> |
call calc_EAM_prepair_rho(i, j, d, r, rijsq ) |
1146 |
> |
endif |
1147 |
> |
|
1148 |
> |
end subroutine do_prepair |
1149 |
> |
|
1150 |
> |
|
1151 |
> |
subroutine do_preforce(nlocal,pot) |
1152 |
> |
integer :: nlocal |
1153 |
> |
real( kind = dp ) :: pot |
1154 |
> |
|
1155 |
> |
if (FF_uses_EAM .and. SIM_uses_EAM) then |
1156 |
> |
call calc_EAM_preforce_Frho(nlocal,pot) |
1157 |
> |
endif |
1158 |
> |
|
1159 |
> |
|
1160 |
> |
end subroutine do_preforce |
1161 |
> |
|
1162 |
> |
|
1163 |
> |
subroutine get_interatomic_vector(q_i, q_j, d, r_sq) |
1164 |
> |
|
1165 |
> |
real (kind = dp), dimension(3) :: q_i |
1166 |
> |
real (kind = dp), dimension(3) :: q_j |
1167 |
> |
real ( kind = dp ), intent(out) :: r_sq |
1168 |
> |
real( kind = dp ) :: d(3), scaled(3) |
1169 |
> |
integer i |
1170 |
> |
|
1171 |
> |
d(1:3) = q_j(1:3) - q_i(1:3) |
1172 |
> |
|
1173 |
> |
! Wrap back into periodic box if necessary |
1174 |
> |
if ( SIM_uses_PBC ) then |
1175 |
> |
|
1176 |
> |
if( .not.boxIsOrthorhombic ) then |
1177 |
> |
! calc the scaled coordinates. |
1178 |
> |
|
1179 |
> |
scaled = matmul(HmatInv, d) |
1180 |
> |
|
1181 |
> |
! wrap the scaled coordinates |
1182 |
> |
|
1183 |
> |
scaled = scaled - anint(scaled) |
1184 |
> |
|
1185 |
> |
|
1186 |
> |
! calc the wrapped real coordinates from the wrapped scaled |
1187 |
> |
! coordinates |
1188 |
> |
|
1189 |
> |
d = matmul(Hmat,scaled) |
1190 |
> |
|
1191 |
> |
else |
1192 |
> |
! calc the scaled coordinates. |
1193 |
> |
|
1194 |
> |
do i = 1, 3 |
1195 |
> |
scaled(i) = d(i) * HmatInv(i,i) |
1196 |
> |
|
1197 |
> |
! wrap the scaled coordinates |
1198 |
> |
|
1199 |
> |
scaled(i) = scaled(i) - anint(scaled(i)) |
1200 |
> |
|
1201 |
> |
! calc the wrapped real coordinates from the wrapped scaled |
1202 |
> |
! coordinates |
1203 |
> |
|
1204 |
> |
d(i) = scaled(i)*Hmat(i,i) |
1205 |
> |
enddo |
1206 |
> |
endif |
1207 |
> |
|
1208 |
> |
endif |
1209 |
> |
|
1210 |
> |
r_sq = dot_product(d,d) |
1211 |
> |
|
1212 |
> |
end subroutine get_interatomic_vector |
1213 |
> |
|
1214 |
> |
subroutine zero_work_arrays() |
1215 |
> |
|
1216 |
|
#ifdef IS_MPI |
1217 |
< |
!! in MPI, we have to look up the unique IDs for each atom |
1218 |
< |
unique_id_1 = AtomRowToGlobal(atom1) |
1219 |
< |
#else |
1220 |
< |
!! in the normal loop, the atom numbers are unique |
1221 |
< |
unique_id_1 = atom1 |
1217 |
> |
|
1218 |
> |
q_Row = 0.0_dp |
1219 |
> |
q_Col = 0.0_dp |
1220 |
> |
|
1221 |
> |
q_group_Row = 0.0_dp |
1222 |
> |
q_group_Col = 0.0_dp |
1223 |
> |
|
1224 |
> |
eFrame_Row = 0.0_dp |
1225 |
> |
eFrame_Col = 0.0_dp |
1226 |
> |
|
1227 |
> |
A_Row = 0.0_dp |
1228 |
> |
A_Col = 0.0_dp |
1229 |
> |
|
1230 |
> |
f_Row = 0.0_dp |
1231 |
> |
f_Col = 0.0_dp |
1232 |
> |
f_Temp = 0.0_dp |
1233 |
> |
|
1234 |
> |
t_Row = 0.0_dp |
1235 |
> |
t_Col = 0.0_dp |
1236 |
> |
t_Temp = 0.0_dp |
1237 |
> |
|
1238 |
> |
pot_Row = 0.0_dp |
1239 |
> |
pot_Col = 0.0_dp |
1240 |
> |
pot_Temp = 0.0_dp |
1241 |
> |
|
1242 |
> |
rf_Row = 0.0_dp |
1243 |
> |
rf_Col = 0.0_dp |
1244 |
> |
rf_Temp = 0.0_dp |
1245 |
> |
|
1246 |
|
#endif |
1247 |
< |
|
1248 |
< |
!! We were called with only one atom, so just check the global exclude |
1249 |
< |
!! list for this atom |
1250 |
< |
if (.not. present(atom2)) then |
1251 |
< |
do i = 1, nExcludes_global |
1252 |
< |
if (excludesGlobal(i) == unique_id_1) then |
1253 |
< |
skip_it = .true. |
1254 |
< |
return |
1255 |
< |
end if |
1256 |
< |
end do |
1257 |
< |
return |
1258 |
< |
end if |
1259 |
< |
|
1247 |
> |
|
1248 |
> |
if (FF_uses_EAM .and. SIM_uses_EAM) then |
1249 |
> |
call clean_EAM() |
1250 |
> |
endif |
1251 |
> |
|
1252 |
> |
rf = 0.0_dp |
1253 |
> |
tau_Temp = 0.0_dp |
1254 |
> |
virial_Temp = 0.0_dp |
1255 |
> |
end subroutine zero_work_arrays |
1256 |
> |
|
1257 |
> |
function skipThisPair(atom1, atom2) result(skip_it) |
1258 |
> |
integer, intent(in) :: atom1 |
1259 |
> |
integer, intent(in), optional :: atom2 |
1260 |
> |
logical :: skip_it |
1261 |
> |
integer :: unique_id_1, unique_id_2 |
1262 |
> |
integer :: me_i,me_j |
1263 |
> |
integer :: i |
1264 |
> |
|
1265 |
> |
skip_it = .false. |
1266 |
> |
|
1267 |
> |
!! there are a number of reasons to skip a pair or a particle |
1268 |
> |
!! mostly we do this to exclude atoms who are involved in short |
1269 |
> |
!! range interactions (bonds, bends, torsions), but we also need |
1270 |
> |
!! to exclude some overcounted interactions that result from |
1271 |
> |
!! the parallel decomposition |
1272 |
> |
|
1273 |
|
#ifdef IS_MPI |
1274 |
< |
unique_id_2 = AtomColToGlobal(atom2) |
1274 |
> |
!! in MPI, we have to look up the unique IDs for each atom |
1275 |
> |
unique_id_1 = AtomRowToGlobal(atom1) |
1276 |
|
#else |
1277 |
< |
unique_id_2 = atom2 |
1277 |
> |
!! in the normal loop, the atom numbers are unique |
1278 |
> |
unique_id_1 = atom1 |
1279 |
|
#endif |
1280 |
< |
|
1280 |
> |
|
1281 |
> |
!! We were called with only one atom, so just check the global exclude |
1282 |
> |
!! list for this atom |
1283 |
> |
if (.not. present(atom2)) then |
1284 |
> |
do i = 1, nExcludes_global |
1285 |
> |
if (excludesGlobal(i) == unique_id_1) then |
1286 |
> |
skip_it = .true. |
1287 |
> |
return |
1288 |
> |
end if |
1289 |
> |
end do |
1290 |
> |
return |
1291 |
> |
end if |
1292 |
> |
|
1293 |
|
#ifdef IS_MPI |
1294 |
< |
!! this situation should only arise in MPI simulations |
1295 |
< |
if (unique_id_1 == unique_id_2) then |
1296 |
< |
skip_it = .true. |
1197 |
< |
return |
1198 |
< |
end if |
1199 |
< |
|
1200 |
< |
!! this prevents us from doing the pair on multiple processors |
1201 |
< |
if (unique_id_1 < unique_id_2) then |
1202 |
< |
if (mod(unique_id_1 + unique_id_2,2) == 0) then |
1203 |
< |
skip_it = .true. |
1204 |
< |
return |
1205 |
< |
endif |
1206 |
< |
else |
1207 |
< |
if (mod(unique_id_1 + unique_id_2,2) == 1) then |
1208 |
< |
skip_it = .true. |
1209 |
< |
return |
1210 |
< |
endif |
1211 |
< |
endif |
1294 |
> |
unique_id_2 = AtomColToGlobal(atom2) |
1295 |
> |
#else |
1296 |
> |
unique_id_2 = atom2 |
1297 |
|
#endif |
1298 |
< |
|
1299 |
< |
!! the rest of these situations can happen in all simulations: |
1300 |
< |
do i = 1, nExcludes_global |
1301 |
< |
if ((excludesGlobal(i) == unique_id_1) .or. & |
1302 |
< |
(excludesGlobal(i) == unique_id_2)) then |
1303 |
< |
skip_it = .true. |
1304 |
< |
return |
1305 |
< |
endif |
1306 |
< |
enddo |
1307 |
< |
|
1308 |
< |
do i = 1, nSkipsForAtom(atom1) |
1309 |
< |
if (skipsForAtom(atom1, i) .eq. unique_id_2) then |
1310 |
< |
skip_it = .true. |
1311 |
< |
return |
1312 |
< |
endif |
1313 |
< |
end do |
1314 |
< |
|
1315 |
< |
return |
1316 |
< |
end function skipThisPair |
1317 |
< |
|
1318 |
< |
function FF_UsesDirectionalAtoms() result(doesit) |
1319 |
< |
logical :: doesit |
1320 |
< |
doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. & |
1321 |
< |
FF_uses_Quadrupoles .or. FF_uses_Sticky .or. & |
1322 |
< |
FF_uses_GayBerne .or. FF_uses_Shapes |
1323 |
< |
end function FF_UsesDirectionalAtoms |
1324 |
< |
|
1325 |
< |
function FF_RequiresPrepairCalc() result(doesit) |
1326 |
< |
logical :: doesit |
1327 |
< |
doesit = FF_uses_EAM |
1328 |
< |
end function FF_RequiresPrepairCalc |
1329 |
< |
|
1330 |
< |
function FF_RequiresPostpairCalc() result(doesit) |
1331 |
< |
logical :: doesit |
1332 |
< |
doesit = FF_uses_RF |
1333 |
< |
end function FF_RequiresPostpairCalc |
1334 |
< |
|
1298 |
> |
|
1299 |
> |
#ifdef IS_MPI |
1300 |
> |
!! this situation should only arise in MPI simulations |
1301 |
> |
if (unique_id_1 == unique_id_2) then |
1302 |
> |
skip_it = .true. |
1303 |
> |
return |
1304 |
> |
end if |
1305 |
> |
|
1306 |
> |
!! this prevents us from doing the pair on multiple processors |
1307 |
> |
if (unique_id_1 < unique_id_2) then |
1308 |
> |
if (mod(unique_id_1 + unique_id_2,2) == 0) then |
1309 |
> |
skip_it = .true. |
1310 |
> |
return |
1311 |
> |
endif |
1312 |
> |
else |
1313 |
> |
if (mod(unique_id_1 + unique_id_2,2) == 1) then |
1314 |
> |
skip_it = .true. |
1315 |
> |
return |
1316 |
> |
endif |
1317 |
> |
endif |
1318 |
> |
#endif |
1319 |
> |
|
1320 |
> |
!! the rest of these situations can happen in all simulations: |
1321 |
> |
do i = 1, nExcludes_global |
1322 |
> |
if ((excludesGlobal(i) == unique_id_1) .or. & |
1323 |
> |
(excludesGlobal(i) == unique_id_2)) then |
1324 |
> |
skip_it = .true. |
1325 |
> |
return |
1326 |
> |
endif |
1327 |
> |
enddo |
1328 |
> |
|
1329 |
> |
do i = 1, nSkipsForAtom(atom1) |
1330 |
> |
if (skipsForAtom(atom1, i) .eq. unique_id_2) then |
1331 |
> |
skip_it = .true. |
1332 |
> |
return |
1333 |
> |
endif |
1334 |
> |
end do |
1335 |
> |
|
1336 |
> |
return |
1337 |
> |
end function skipThisPair |
1338 |
> |
|
1339 |
> |
function FF_UsesDirectionalAtoms() result(doesit) |
1340 |
> |
logical :: doesit |
1341 |
> |
doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. & |
1342 |
> |
FF_uses_Quadrupoles .or. FF_uses_Sticky .or. & |
1343 |
> |
FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes |
1344 |
> |
end function FF_UsesDirectionalAtoms |
1345 |
> |
|
1346 |
> |
function FF_RequiresPrepairCalc() result(doesit) |
1347 |
> |
logical :: doesit |
1348 |
> |
doesit = FF_uses_EAM |
1349 |
> |
end function FF_RequiresPrepairCalc |
1350 |
> |
|
1351 |
> |
function FF_RequiresPostpairCalc() result(doesit) |
1352 |
> |
logical :: doesit |
1353 |
> |
doesit = FF_uses_RF |
1354 |
> |
end function FF_RequiresPostpairCalc |
1355 |
> |
|
1356 |
|
#ifdef PROFILE |
1357 |
< |
function getforcetime() result(totalforcetime) |
1358 |
< |
real(kind=dp) :: totalforcetime |
1359 |
< |
totalforcetime = forcetime |
1360 |
< |
end function getforcetime |
1357 |
> |
function getforcetime() result(totalforcetime) |
1358 |
> |
real(kind=dp) :: totalforcetime |
1359 |
> |
totalforcetime = forcetime |
1360 |
> |
end function getforcetime |
1361 |
|
#endif |
1256 |
– |
|
1257 |
– |
!! This cleans componets of force arrays belonging only to fortran |
1362 |
|
|
1363 |
< |
subroutine add_stress_tensor(dpair, fpair) |
1364 |
< |
|
1365 |
< |
real( kind = dp ), dimension(3), intent(in) :: dpair, fpair |
1366 |
< |
|
1367 |
< |
! because the d vector is the rj - ri vector, and |
1368 |
< |
! because fx, fy, fz are the force on atom i, we need a |
1369 |
< |
! negative sign here: |
1370 |
< |
|
1371 |
< |
tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1) |
1372 |
< |
tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2) |
1373 |
< |
tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3) |
1374 |
< |
tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1) |
1375 |
< |
tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2) |
1376 |
< |
tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3) |
1377 |
< |
tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1) |
1378 |
< |
tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2) |
1379 |
< |
tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3) |
1380 |
< |
|
1381 |
< |
virial_Temp = virial_Temp + & |
1382 |
< |
(tau_Temp(1) + tau_Temp(5) + tau_Temp(9)) |
1383 |
< |
|
1384 |
< |
end subroutine add_stress_tensor |
1385 |
< |
|
1363 |
> |
!! This cleans componets of force arrays belonging only to fortran |
1364 |
> |
|
1365 |
> |
subroutine add_stress_tensor(dpair, fpair) |
1366 |
> |
|
1367 |
> |
real( kind = dp ), dimension(3), intent(in) :: dpair, fpair |
1368 |
> |
|
1369 |
> |
! because the d vector is the rj - ri vector, and |
1370 |
> |
! because fx, fy, fz are the force on atom i, we need a |
1371 |
> |
! negative sign here: |
1372 |
> |
|
1373 |
> |
tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1) |
1374 |
> |
tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2) |
1375 |
> |
tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3) |
1376 |
> |
tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1) |
1377 |
> |
tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2) |
1378 |
> |
tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3) |
1379 |
> |
tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1) |
1380 |
> |
tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2) |
1381 |
> |
tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3) |
1382 |
> |
|
1383 |
> |
virial_Temp = virial_Temp + & |
1384 |
> |
(tau_Temp(1) + tau_Temp(5) + tau_Temp(9)) |
1385 |
> |
|
1386 |
> |
end subroutine add_stress_tensor |
1387 |
> |
|
1388 |
|
end module doForces |