45 |
|
|
46 |
|
!! @author Charles F. Vardeman II |
47 |
|
!! @author Matthew Meineke |
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 $ |
48 |
> |
!! @version $Id: doForces.F90,v 1.69 2005-11-21 22:58:35 gezelter Exp $, $Date: 2005-11-21 22:58:35 $, $Name: not supported by cvs2svn $, $Revision: 1.69 $ |
49 |
|
|
50 |
|
|
51 |
|
module doForces |
58 |
|
use lj |
59 |
|
use sticky |
60 |
|
use electrostatic_module |
61 |
< |
use reaction_field |
62 |
< |
use gb_pair |
61 |
> |
use gayberne |
62 |
|
use shapes |
63 |
|
use vector_class |
64 |
|
use eam |
65 |
+ |
use suttonchen |
66 |
|
use status |
67 |
|
#ifdef IS_MPI |
68 |
|
use mpiSimulation |
73 |
|
|
74 |
|
#define __FORTRAN90 |
75 |
|
#include "UseTheForce/fSwitchingFunction.h" |
76 |
+ |
#include "UseTheForce/fCutoffPolicy.h" |
77 |
|
#include "UseTheForce/DarkSide/fInteractionMap.h" |
78 |
+ |
#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" |
79 |
|
|
80 |
+ |
|
81 |
|
INTEGER, PARAMETER:: PREPAIR_LOOP = 1 |
82 |
|
INTEGER, PARAMETER:: PAIR_LOOP = 2 |
83 |
|
|
81 |
– |
logical, save :: haveRlist = .false. |
84 |
|
logical, save :: haveNeighborList = .false. |
85 |
|
logical, save :: haveSIMvariables = .false. |
86 |
|
logical, save :: haveSaneForceField = .false. |
87 |
< |
logical, save :: haveInteractionMap = .false. |
87 |
> |
logical, save :: haveInteractionHash = .false. |
88 |
> |
logical, save :: haveGtypeCutoffMap = .false. |
89 |
> |
logical, save :: haveDefaultCutoffs = .false. |
90 |
> |
logical, save :: haveSkinThickness = .false. |
91 |
> |
logical, save :: haveElectrostaticSummationMethod = .false. |
92 |
> |
logical, save :: haveCutoffPolicy = .false. |
93 |
> |
logical, save :: VisitCutoffsAfterComputing = .false. |
94 |
|
|
95 |
|
logical, save :: FF_uses_DirectionalAtoms |
88 |
– |
logical, save :: FF_uses_LennardJones |
89 |
– |
logical, save :: FF_uses_Electrostatics |
90 |
– |
logical, save :: FF_uses_Charges |
96 |
|
logical, save :: FF_uses_Dipoles |
92 |
– |
logical, save :: FF_uses_Quadrupoles |
93 |
– |
logical, save :: FF_uses_Sticky |
94 |
– |
logical, save :: FF_uses_StickyPower |
97 |
|
logical, save :: FF_uses_GayBerne |
98 |
|
logical, save :: FF_uses_EAM |
99 |
< |
logical, save :: FF_uses_Shapes |
100 |
< |
logical, save :: FF_uses_FLARB |
101 |
< |
logical, save :: FF_uses_RF |
99 |
> |
logical, save :: FF_uses_SC |
100 |
> |
logical, save :: FF_uses_MEAM |
101 |
> |
|
102 |
|
|
103 |
|
logical, save :: SIM_uses_DirectionalAtoms |
102 |
– |
logical, save :: SIM_uses_LennardJones |
103 |
– |
logical, save :: SIM_uses_Electrostatics |
104 |
– |
logical, save :: SIM_uses_Charges |
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 |
104 |
|
logical, save :: SIM_uses_EAM |
105 |
< |
logical, save :: SIM_uses_Shapes |
106 |
< |
logical, save :: SIM_uses_FLARB |
113 |
< |
logical, save :: SIM_uses_RF |
105 |
> |
logical, save :: SIM_uses_SC |
106 |
> |
logical, save :: SIM_uses_MEAM |
107 |
|
logical, save :: SIM_requires_postpair_calc |
108 |
|
logical, save :: SIM_requires_prepair_calc |
109 |
|
logical, save :: SIM_uses_PBC |
117 |
– |
logical, save :: SIM_uses_molecular_cutoffs |
110 |
|
|
111 |
< |
!!!GO AWAY--------- |
112 |
< |
!!!!!real(kind=dp), save :: rlist, rlistsq |
111 |
> |
integer, save :: electrostaticSummationMethod |
112 |
> |
integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY |
113 |
|
|
114 |
+ |
real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut |
115 |
+ |
real(kind=dp), save :: skinThickness |
116 |
+ |
logical, save :: defaultDoShift |
117 |
+ |
|
118 |
|
public :: init_FF |
119 |
+ |
public :: setCutoffs |
120 |
+ |
public :: cWasLame |
121 |
+ |
public :: setElectrostaticMethod |
122 |
+ |
public :: setCutoffPolicy |
123 |
+ |
public :: setSkinThickness |
124 |
|
public :: do_force_loop |
124 |
– |
! public :: setRlistDF |
125 |
– |
!public :: addInteraction |
126 |
– |
!public :: setInteractionHash |
127 |
– |
!public :: getInteractionHash |
128 |
– |
public :: createInteractionMap |
129 |
– |
public :: createRcuts |
125 |
|
|
126 |
|
#ifdef PROFILE |
127 |
|
public :: getforcetime |
129 |
|
real :: forceTimeInitial, forceTimeFinal |
130 |
|
integer :: nLoops |
131 |
|
#endif |
137 |
– |
|
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 |
132 |
|
|
133 |
< |
type(Interaction), dimension(:,:),allocatable :: InteractionMap |
134 |
< |
|
133 |
> |
!! Variables for cutoff mapping and interaction mapping |
134 |
> |
! Bit hash to determine pair-pair interactions. |
135 |
> |
integer, dimension(:,:), allocatable :: InteractionHash |
136 |
> |
real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff |
137 |
> |
real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow |
138 |
> |
real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol |
139 |
|
|
140 |
< |
|
140 |
> |
integer, dimension(:), allocatable, target :: groupToGtypeRow |
141 |
> |
integer, dimension(:), pointer :: groupToGtypeCol => null() |
142 |
> |
|
143 |
> |
real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow |
144 |
> |
real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol |
145 |
> |
type ::gtypeCutoffs |
146 |
> |
real(kind=dp) :: rcut |
147 |
> |
real(kind=dp) :: rcutsq |
148 |
> |
real(kind=dp) :: rlistsq |
149 |
> |
end type gtypeCutoffs |
150 |
> |
type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap |
151 |
> |
|
152 |
|
contains |
153 |
|
|
154 |
< |
|
151 |
< |
subroutine createInteractionMap(status) |
154 |
> |
subroutine createInteractionHash() |
155 |
|
integer :: nAtypes |
153 |
– |
integer :: status |
156 |
|
integer :: i |
157 |
|
integer :: j |
158 |
< |
integer :: ihash |
159 |
< |
real(kind=dp) :: myRcut |
158 |
< |
! Test Types |
158 |
> |
integer :: iHash |
159 |
> |
!! Test Types |
160 |
|
logical :: i_is_LJ |
161 |
|
logical :: i_is_Elect |
162 |
|
logical :: i_is_Sticky |
164 |
|
logical :: i_is_GB |
165 |
|
logical :: i_is_EAM |
166 |
|
logical :: i_is_Shape |
167 |
+ |
logical :: i_is_SC |
168 |
+ |
logical :: i_is_MEAM |
169 |
|
logical :: j_is_LJ |
170 |
|
logical :: j_is_Elect |
171 |
|
logical :: j_is_Sticky |
173 |
|
logical :: j_is_GB |
174 |
|
logical :: j_is_EAM |
175 |
|
logical :: j_is_Shape |
176 |
< |
|
177 |
< |
|
176 |
> |
logical :: j_is_SC |
177 |
> |
logical :: j_is_MEAM |
178 |
> |
real(kind=dp) :: myRcut |
179 |
> |
|
180 |
|
if (.not. associated(atypes)) then |
181 |
< |
call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!") |
177 |
< |
status = -1 |
181 |
> |
call handleError("doForces", "atypes was not present before call of createInteractionHash!") |
182 |
|
return |
183 |
|
endif |
184 |
|
|
185 |
|
nAtypes = getSize(atypes) |
186 |
|
|
187 |
|
if (nAtypes == 0) then |
188 |
< |
status = -1 |
188 |
> |
call handleError("doForces", "nAtypes was zero during call of createInteractionHash!") |
189 |
|
return |
190 |
|
end if |
191 |
|
|
192 |
< |
if (.not. allocated(InteractionMap)) then |
193 |
< |
allocate(InteractionMap(nAtypes,nAtypes)) |
192 |
> |
if (.not. allocated(InteractionHash)) then |
193 |
> |
allocate(InteractionHash(nAtypes,nAtypes)) |
194 |
> |
else |
195 |
> |
deallocate(InteractionHash) |
196 |
> |
allocate(InteractionHash(nAtypes,nAtypes)) |
197 |
|
endif |
198 |
+ |
|
199 |
+ |
if (.not. allocated(atypeMaxCutoff)) then |
200 |
+ |
allocate(atypeMaxCutoff(nAtypes)) |
201 |
+ |
else |
202 |
+ |
deallocate(atypeMaxCutoff) |
203 |
+ |
allocate(atypeMaxCutoff(nAtypes)) |
204 |
+ |
endif |
205 |
|
|
206 |
|
do i = 1, nAtypes |
207 |
|
call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) |
211 |
|
call getElementProperty(atypes, i, "is_GayBerne", i_is_GB) |
212 |
|
call getElementProperty(atypes, i, "is_EAM", i_is_EAM) |
213 |
|
call getElementProperty(atypes, i, "is_Shape", i_is_Shape) |
214 |
+ |
call getElementProperty(atypes, i, "is_SC", i_is_SC) |
215 |
+ |
call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM) |
216 |
|
|
217 |
|
do j = i, nAtypes |
218 |
|
|
226 |
|
call getElementProperty(atypes, j, "is_GayBerne", j_is_GB) |
227 |
|
call getElementProperty(atypes, j, "is_EAM", j_is_EAM) |
228 |
|
call getElementProperty(atypes, j, "is_Shape", j_is_Shape) |
229 |
+ |
call getElementProperty(atypes, j, "is_SC", j_is_SC) |
230 |
+ |
call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM) |
231 |
|
|
232 |
|
if (i_is_LJ .and. j_is_LJ) then |
233 |
|
iHash = ior(iHash, LJ_PAIR) |
249 |
|
iHash = ior(iHash, EAM_PAIR) |
250 |
|
endif |
251 |
|
|
252 |
+ |
if (i_is_SC .and. j_is_SC) then |
253 |
+ |
iHash = ior(iHash, SC_PAIR) |
254 |
+ |
endif |
255 |
+ |
|
256 |
|
if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR) |
257 |
|
if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ) |
258 |
|
if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ) |
262 |
|
if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ) |
263 |
|
|
264 |
|
|
265 |
< |
InteractionMap(i,j)%InteractionHash = iHash |
266 |
< |
InteractionMap(j,i)%InteractionHash = iHash |
265 |
> |
InteractionHash(i,j) = iHash |
266 |
> |
InteractionHash(j,i) = iHash |
267 |
|
|
268 |
|
end do |
269 |
|
|
270 |
|
end do |
249 |
– |
end subroutine createInteractionMap |
271 |
|
|
272 |
< |
! 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. |
273 |
< |
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 |
272 |
> |
haveInteractionHash = .true. |
273 |
> |
end subroutine createInteractionHash |
274 |
|
|
275 |
< |
if(.not. allocated(InteractionMap)) return |
275 |
> |
subroutine createGtypeCutoffMap() |
276 |
> |
|
277 |
> |
logical :: i_is_LJ |
278 |
> |
logical :: i_is_Elect |
279 |
> |
logical :: i_is_Sticky |
280 |
> |
logical :: i_is_StickyP |
281 |
> |
logical :: i_is_GB |
282 |
> |
logical :: i_is_EAM |
283 |
> |
logical :: i_is_Shape |
284 |
> |
logical :: GtypeFound |
285 |
> |
|
286 |
> |
integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend |
287 |
> |
integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j |
288 |
> |
integer :: nGroupsInRow |
289 |
> |
integer :: nGroupsInCol |
290 |
> |
integer :: nGroupTypesRow,nGroupTypesCol |
291 |
> |
real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol |
292 |
> |
real(kind=dp) :: biggestAtypeCutoff |
293 |
|
|
294 |
+ |
if (.not. haveInteractionHash) then |
295 |
+ |
call createInteractionHash() |
296 |
+ |
endif |
297 |
+ |
#ifdef IS_MPI |
298 |
+ |
nGroupsInRow = getNgroupsInRow(plan_group_row) |
299 |
+ |
nGroupsInCol = getNgroupsInCol(plan_group_col) |
300 |
+ |
#endif |
301 |
|
nAtypes = getSize(atypes) |
302 |
< |
! If we pass a default rcut, set all atypes to that cutoff distance |
303 |
< |
if(present(defaultRList)) then |
304 |
< |
InteractionMap(:,:)%rList = defaultRList |
305 |
< |
InteractionMap(:,:)%rListSq = defaultRList*defaultRList |
306 |
< |
haveRlist = .true. |
307 |
< |
return |
308 |
< |
end if |
309 |
< |
|
310 |
< |
do map_i = 1,nAtypes |
311 |
< |
do map_j = map_i,nAtypes |
312 |
< |
iMap = InteractionMap(map_i, map_j)%InteractionHash |
302 |
> |
! Set all of the initial cutoffs to zero. |
303 |
> |
atypeMaxCutoff = 0.0_dp |
304 |
> |
do i = 1, nAtypes |
305 |
> |
if (SimHasAtype(i)) then |
306 |
> |
call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) |
307 |
> |
call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect) |
308 |
> |
call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky) |
309 |
> |
call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP) |
310 |
> |
call getElementProperty(atypes, i, "is_GayBerne", i_is_GB) |
311 |
> |
call getElementProperty(atypes, i, "is_EAM", i_is_EAM) |
312 |
> |
call getElementProperty(atypes, i, "is_Shape", i_is_Shape) |
313 |
|
|
314 |
< |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
315 |
< |
! thisRCut = getLJCutOff(map_i,map_j) |
316 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
314 |
> |
|
315 |
> |
if (haveDefaultCutoffs) then |
316 |
> |
atypeMaxCutoff(i) = defaultRcut |
317 |
> |
else |
318 |
> |
if (i_is_LJ) then |
319 |
> |
thisRcut = getSigma(i) * 2.5_dp |
320 |
> |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
321 |
> |
endif |
322 |
> |
if (i_is_Elect) then |
323 |
> |
thisRcut = defaultRcut |
324 |
> |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
325 |
> |
endif |
326 |
> |
if (i_is_Sticky) then |
327 |
> |
thisRcut = getStickyCut(i) |
328 |
> |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
329 |
> |
endif |
330 |
> |
if (i_is_StickyP) then |
331 |
> |
thisRcut = getStickyPowerCut(i) |
332 |
> |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
333 |
> |
endif |
334 |
> |
if (i_is_GB) then |
335 |
> |
thisRcut = getGayBerneCut(i) |
336 |
> |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
337 |
> |
endif |
338 |
> |
if (i_is_EAM) then |
339 |
> |
thisRcut = getEAMCut(i) |
340 |
> |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
341 |
> |
endif |
342 |
> |
if (i_is_Shape) then |
343 |
> |
thisRcut = getShapeCut(i) |
344 |
> |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
345 |
> |
endif |
346 |
|
endif |
347 |
< |
|
348 |
< |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
349 |
< |
! thisRCut = getElectrostaticCutOff(map_i,map_j) |
282 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
347 |
> |
|
348 |
> |
if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then |
349 |
> |
biggestAtypeCutoff = atypeMaxCutoff(i) |
350 |
|
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 |
351 |
|
|
352 |
+ |
endif |
353 |
+ |
enddo |
354 |
+ |
|
355 |
+ |
istart = 1 |
356 |
+ |
jstart = 1 |
357 |
+ |
#ifdef IS_MPI |
358 |
+ |
iend = nGroupsInRow |
359 |
+ |
jend = nGroupsInCol |
360 |
+ |
#else |
361 |
+ |
iend = nGroups |
362 |
+ |
jend = nGroups |
363 |
+ |
#endif |
364 |
+ |
|
365 |
+ |
!! allocate the groupToGtype and gtypeMaxCutoff here. |
366 |
+ |
if(.not.allocated(groupToGtypeRow)) then |
367 |
+ |
! allocate(groupToGtype(iend)) |
368 |
+ |
allocate(groupToGtypeRow(iend)) |
369 |
+ |
else |
370 |
+ |
deallocate(groupToGtypeRow) |
371 |
+ |
allocate(groupToGtypeRow(iend)) |
372 |
+ |
endif |
373 |
+ |
if(.not.allocated(groupMaxCutoffRow)) then |
374 |
+ |
allocate(groupMaxCutoffRow(iend)) |
375 |
+ |
else |
376 |
+ |
deallocate(groupMaxCutoffRow) |
377 |
+ |
allocate(groupMaxCutoffRow(iend)) |
378 |
+ |
end if |
379 |
|
|
380 |
< |
!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF |
381 |
< |
!!$ subroutine setRlistDF( this_rlist ) |
382 |
< |
!!$ |
383 |
< |
!!$ real(kind=dp) :: this_rlist |
384 |
< |
!!$ |
385 |
< |
!!$ rlist = this_rlist |
333 |
< |
!!$ rlistsq = rlist * rlist |
334 |
< |
!!$ |
335 |
< |
!!$ haveRlist = .true. |
336 |
< |
!!$ |
337 |
< |
!!$ end subroutine setRlistDF |
380 |
> |
if(.not.allocated(gtypeMaxCutoffRow)) then |
381 |
> |
allocate(gtypeMaxCutoffRow(iend)) |
382 |
> |
else |
383 |
> |
deallocate(gtypeMaxCutoffRow) |
384 |
> |
allocate(gtypeMaxCutoffRow(iend)) |
385 |
> |
endif |
386 |
|
|
387 |
|
|
388 |
< |
subroutine setSimVariables() |
389 |
< |
SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms() |
390 |
< |
SIM_uses_LennardJones = SimUsesLennardJones() |
391 |
< |
SIM_uses_Electrostatics = SimUsesElectrostatics() |
392 |
< |
SIM_uses_Charges = SimUsesCharges() |
393 |
< |
SIM_uses_Dipoles = SimUsesDipoles() |
394 |
< |
SIM_uses_Sticky = SimUsesSticky() |
395 |
< |
SIM_uses_StickyPower = SimUsesStickyPower() |
348 |
< |
SIM_uses_GayBerne = SimUsesGayBerne() |
349 |
< |
SIM_uses_EAM = SimUsesEAM() |
350 |
< |
SIM_uses_Shapes = SimUsesShapes() |
351 |
< |
SIM_uses_FLARB = SimUsesFLARB() |
352 |
< |
SIM_uses_RF = SimUsesRF() |
353 |
< |
SIM_requires_postpair_calc = SimRequiresPostpairCalc() |
354 |
< |
SIM_requires_prepair_calc = SimRequiresPrepairCalc() |
355 |
< |
SIM_uses_PBC = SimUsesPBC() |
388 |
> |
#ifdef IS_MPI |
389 |
> |
! We only allocate new storage if we are in MPI because Ncol /= Nrow |
390 |
> |
if(.not.associated(groupToGtypeCol)) then |
391 |
> |
allocate(groupToGtypeCol(jend)) |
392 |
> |
else |
393 |
> |
deallocate(groupToGtypeCol) |
394 |
> |
allocate(groupToGtypeCol(jend)) |
395 |
> |
end if |
396 |
|
|
397 |
< |
haveSIMvariables = .true. |
397 |
> |
if(.not.associated(groupToGtypeCol)) then |
398 |
> |
allocate(groupToGtypeCol(jend)) |
399 |
> |
else |
400 |
> |
deallocate(groupToGtypeCol) |
401 |
> |
allocate(groupToGtypeCol(jend)) |
402 |
> |
end if |
403 |
> |
if(.not.associated(gtypeMaxCutoffCol)) then |
404 |
> |
allocate(gtypeMaxCutoffCol(jend)) |
405 |
> |
else |
406 |
> |
deallocate(gtypeMaxCutoffCol) |
407 |
> |
allocate(gtypeMaxCutoffCol(jend)) |
408 |
> |
end if |
409 |
|
|
410 |
< |
return |
411 |
< |
end subroutine setSimVariables |
410 |
> |
groupMaxCutoffCol = 0.0_dp |
411 |
> |
gtypeMaxCutoffCol = 0.0_dp |
412 |
|
|
413 |
+ |
#endif |
414 |
+ |
groupMaxCutoffRow = 0.0_dp |
415 |
+ |
gtypeMaxCutoffRow = 0.0_dp |
416 |
+ |
|
417 |
+ |
|
418 |
+ |
!! first we do a single loop over the cutoff groups to find the |
419 |
+ |
!! largest cutoff for any atypes present in this group. We also |
420 |
+ |
!! create gtypes at this point. |
421 |
+ |
|
422 |
+ |
tol = 1.0d-6 |
423 |
+ |
nGroupTypesRow = 0 |
424 |
+ |
|
425 |
+ |
do i = istart, iend |
426 |
+ |
n_in_i = groupStartRow(i+1) - groupStartRow(i) |
427 |
+ |
groupMaxCutoffRow(i) = 0.0_dp |
428 |
+ |
do ia = groupStartRow(i), groupStartRow(i+1)-1 |
429 |
+ |
atom1 = groupListRow(ia) |
430 |
+ |
#ifdef IS_MPI |
431 |
+ |
me_i = atid_row(atom1) |
432 |
+ |
#else |
433 |
+ |
me_i = atid(atom1) |
434 |
+ |
#endif |
435 |
+ |
if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then |
436 |
+ |
groupMaxCutoffRow(i)=atypeMaxCutoff(me_i) |
437 |
+ |
endif |
438 |
+ |
enddo |
439 |
+ |
|
440 |
+ |
if (nGroupTypesRow.eq.0) then |
441 |
+ |
nGroupTypesRow = nGroupTypesRow + 1 |
442 |
+ |
gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i) |
443 |
+ |
groupToGtypeRow(i) = nGroupTypesRow |
444 |
+ |
else |
445 |
+ |
GtypeFound = .false. |
446 |
+ |
do g = 1, nGroupTypesRow |
447 |
+ |
if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then |
448 |
+ |
groupToGtypeRow(i) = g |
449 |
+ |
GtypeFound = .true. |
450 |
+ |
endif |
451 |
+ |
enddo |
452 |
+ |
if (.not.GtypeFound) then |
453 |
+ |
nGroupTypesRow = nGroupTypesRow + 1 |
454 |
+ |
gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i) |
455 |
+ |
groupToGtypeRow(i) = nGroupTypesRow |
456 |
+ |
endif |
457 |
+ |
endif |
458 |
+ |
enddo |
459 |
+ |
|
460 |
+ |
#ifdef IS_MPI |
461 |
+ |
do j = jstart, jend |
462 |
+ |
n_in_j = groupStartCol(j+1) - groupStartCol(j) |
463 |
+ |
groupMaxCutoffCol(j) = 0.0_dp |
464 |
+ |
do ja = groupStartCol(j), groupStartCol(j+1)-1 |
465 |
+ |
atom1 = groupListCol(ja) |
466 |
+ |
|
467 |
+ |
me_j = atid_col(atom1) |
468 |
+ |
|
469 |
+ |
if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then |
470 |
+ |
groupMaxCutoffCol(j)=atypeMaxCutoff(me_j) |
471 |
+ |
endif |
472 |
+ |
enddo |
473 |
+ |
|
474 |
+ |
if (nGroupTypesCol.eq.0) then |
475 |
+ |
nGroupTypesCol = nGroupTypesCol + 1 |
476 |
+ |
gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j) |
477 |
+ |
groupToGtypeCol(j) = nGroupTypesCol |
478 |
+ |
else |
479 |
+ |
GtypeFound = .false. |
480 |
+ |
do g = 1, nGroupTypesCol |
481 |
+ |
if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then |
482 |
+ |
groupToGtypeCol(j) = g |
483 |
+ |
GtypeFound = .true. |
484 |
+ |
endif |
485 |
+ |
enddo |
486 |
+ |
if (.not.GtypeFound) then |
487 |
+ |
nGroupTypesCol = nGroupTypesCol + 1 |
488 |
+ |
gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j) |
489 |
+ |
groupToGtypeCol(j) = nGroupTypesCol |
490 |
+ |
endif |
491 |
+ |
endif |
492 |
+ |
enddo |
493 |
+ |
|
494 |
+ |
#else |
495 |
+ |
! Set pointers to information we just found |
496 |
+ |
nGroupTypesCol = nGroupTypesRow |
497 |
+ |
groupToGtypeCol => groupToGtypeRow |
498 |
+ |
gtypeMaxCutoffCol => gtypeMaxCutoffRow |
499 |
+ |
groupMaxCutoffCol => groupMaxCutoffRow |
500 |
+ |
#endif |
501 |
+ |
|
502 |
+ |
!! allocate the gtypeCutoffMap here. |
503 |
+ |
allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol)) |
504 |
+ |
!! then we do a double loop over all the group TYPES to find the cutoff |
505 |
+ |
!! map between groups of two types |
506 |
+ |
tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol)) |
507 |
+ |
|
508 |
+ |
do i = 1, nGroupTypesRow |
509 |
+ |
do j = 1, nGroupTypesCol |
510 |
+ |
|
511 |
+ |
select case(cutoffPolicy) |
512 |
+ |
case(TRADITIONAL_CUTOFF_POLICY) |
513 |
+ |
thisRcut = tradRcut |
514 |
+ |
case(MIX_CUTOFF_POLICY) |
515 |
+ |
thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j)) |
516 |
+ |
case(MAX_CUTOFF_POLICY) |
517 |
+ |
thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j)) |
518 |
+ |
case default |
519 |
+ |
call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy") |
520 |
+ |
return |
521 |
+ |
end select |
522 |
+ |
gtypeCutoffMap(i,j)%rcut = thisRcut |
523 |
+ |
|
524 |
+ |
if (thisRcut.gt.largestRcut) largestRcut = thisRcut |
525 |
+ |
|
526 |
+ |
gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut |
527 |
+ |
|
528 |
+ |
if (.not.haveSkinThickness) then |
529 |
+ |
skinThickness = 1.0_dp |
530 |
+ |
endif |
531 |
+ |
|
532 |
+ |
gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2 |
533 |
+ |
|
534 |
+ |
! sanity check |
535 |
+ |
|
536 |
+ |
if (haveDefaultCutoffs) then |
537 |
+ |
if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then |
538 |
+ |
call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff") |
539 |
+ |
endif |
540 |
+ |
endif |
541 |
+ |
enddo |
542 |
+ |
enddo |
543 |
+ |
|
544 |
+ |
if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow) |
545 |
+ |
if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow) |
546 |
+ |
if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff) |
547 |
+ |
#ifdef IS_MPI |
548 |
+ |
if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol) |
549 |
+ |
if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol) |
550 |
+ |
#endif |
551 |
+ |
groupMaxCutoffCol => null() |
552 |
+ |
gtypeMaxCutoffCol => null() |
553 |
+ |
|
554 |
+ |
haveGtypeCutoffMap = .true. |
555 |
+ |
end subroutine createGtypeCutoffMap |
556 |
+ |
|
557 |
+ |
subroutine setCutoffs(defRcut, defRsw) |
558 |
+ |
|
559 |
+ |
real(kind=dp),intent(in) :: defRcut, defRsw |
560 |
+ |
character(len = statusMsgSize) :: errMsg |
561 |
+ |
integer :: localError |
562 |
+ |
|
563 |
+ |
defaultRcut = defRcut |
564 |
+ |
defaultRsw = defRsw |
565 |
+ |
|
566 |
+ |
defaultDoShift = .false. |
567 |
+ |
if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then |
568 |
+ |
|
569 |
+ |
write(errMsg, *) & |
570 |
+ |
'cutoffRadius and switchingRadius are set to the same', newline & |
571 |
+ |
// tab, 'value. OOPSE will use shifted ', newline & |
572 |
+ |
// tab, 'potentials instead of switching functions.' |
573 |
+ |
|
574 |
+ |
call handleInfo("setCutoffs", errMsg) |
575 |
+ |
|
576 |
+ |
defaultDoShift = .true. |
577 |
+ |
|
578 |
+ |
endif |
579 |
+ |
|
580 |
+ |
localError = 0 |
581 |
+ |
call setLJDefaultCutoff( defaultRcut, defaultDoShift ) |
582 |
+ |
call setCutoffEAM( defaultRcut, localError) |
583 |
+ |
if (localError /= 0) then |
584 |
+ |
write(errMsg, *) 'An error has occured in setting the EAM cutoff' |
585 |
+ |
call handleError("setCutoffs", errMsg) |
586 |
+ |
end if |
587 |
+ |
call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut) |
588 |
+ |
|
589 |
+ |
haveDefaultCutoffs = .true. |
590 |
+ |
end subroutine setCutoffs |
591 |
+ |
|
592 |
+ |
subroutine cWasLame() |
593 |
+ |
|
594 |
+ |
VisitCutoffsAfterComputing = .true. |
595 |
+ |
return |
596 |
+ |
|
597 |
+ |
end subroutine cWasLame |
598 |
+ |
|
599 |
+ |
subroutine setCutoffPolicy(cutPolicy) |
600 |
+ |
|
601 |
+ |
integer, intent(in) :: cutPolicy |
602 |
+ |
|
603 |
+ |
cutoffPolicy = cutPolicy |
604 |
+ |
haveCutoffPolicy = .true. |
605 |
+ |
|
606 |
+ |
call createGtypeCutoffMap() |
607 |
+ |
|
608 |
+ |
end subroutine setCutoffPolicy |
609 |
+ |
|
610 |
+ |
subroutine setElectrostaticMethod( thisESM ) |
611 |
+ |
|
612 |
+ |
integer, intent(in) :: thisESM |
613 |
+ |
|
614 |
+ |
electrostaticSummationMethod = thisESM |
615 |
+ |
haveElectrostaticSummationMethod = .true. |
616 |
+ |
|
617 |
+ |
end subroutine setElectrostaticMethod |
618 |
+ |
|
619 |
+ |
subroutine setSkinThickness( thisSkin ) |
620 |
+ |
|
621 |
+ |
real(kind=dp), intent(in) :: thisSkin |
622 |
+ |
|
623 |
+ |
skinThickness = thisSkin |
624 |
+ |
haveSkinThickness = .true. |
625 |
+ |
|
626 |
+ |
call createGtypeCutoffMap() |
627 |
+ |
|
628 |
+ |
end subroutine setSkinThickness |
629 |
+ |
|
630 |
+ |
subroutine setSimVariables() |
631 |
+ |
SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms() |
632 |
+ |
SIM_uses_EAM = SimUsesEAM() |
633 |
+ |
SIM_requires_postpair_calc = SimRequiresPostpairCalc() |
634 |
+ |
SIM_requires_prepair_calc = SimRequiresPrepairCalc() |
635 |
+ |
SIM_uses_PBC = SimUsesPBC() |
636 |
+ |
|
637 |
+ |
haveSIMvariables = .true. |
638 |
+ |
|
639 |
+ |
return |
640 |
+ |
end subroutine setSimVariables |
641 |
+ |
|
642 |
|
subroutine doReadyCheck(error) |
643 |
|
integer, intent(out) :: error |
644 |
|
|
646 |
|
|
647 |
|
error = 0 |
648 |
|
|
649 |
< |
if (.not. haveInteractionMap) then |
649 |
> |
if (.not. haveInteractionHash) then |
650 |
> |
call createInteractionHash() |
651 |
> |
endif |
652 |
|
|
653 |
< |
myStatus = 0 |
653 |
> |
if (.not. haveGtypeCutoffMap) then |
654 |
> |
call createGtypeCutoffMap() |
655 |
> |
endif |
656 |
|
|
373 |
– |
call createInteractionMap(myStatus) |
657 |
|
|
658 |
< |
if (myStatus .ne. 0) then |
659 |
< |
write(default_error, *) 'createInteractionMap failed in doForces!' |
377 |
< |
error = -1 |
378 |
< |
return |
379 |
< |
endif |
658 |
> |
if (VisitCutoffsAfterComputing) then |
659 |
> |
call set_switch(GROUP_SWITCH, largestRcut, largestRcut) |
660 |
|
endif |
661 |
|
|
662 |
+ |
|
663 |
|
if (.not. haveSIMvariables) then |
664 |
|
call setSimVariables() |
665 |
|
endif |
666 |
|
|
667 |
< |
if (.not. haveRlist) then |
668 |
< |
write(default_error, *) 'rList has not been set in doForces!' |
669 |
< |
error = -1 |
670 |
< |
return |
671 |
< |
endif |
667 |
> |
! if (.not. haveRlist) then |
668 |
> |
! write(default_error, *) 'rList has not been set in doForces!' |
669 |
> |
! error = -1 |
670 |
> |
! return |
671 |
> |
! endif |
672 |
|
|
673 |
|
if (.not. haveNeighborList) then |
674 |
|
write(default_error, *) 'neighbor list has not been initialized in doForces!' |
693 |
|
end subroutine doReadyCheck |
694 |
|
|
695 |
|
|
696 |
< |
subroutine init_FF(use_RF_c, thisStat) |
696 |
> |
subroutine init_FF(thisStat) |
697 |
|
|
417 |
– |
logical, intent(in) :: use_RF_c |
418 |
– |
|
698 |
|
integer, intent(out) :: thisStat |
699 |
|
integer :: my_status, nMatches |
700 |
|
integer, pointer :: MatchList(:) => null() |
422 |
– |
real(kind=dp) :: rcut, rrf, rt, dielect |
701 |
|
|
702 |
|
!! assume things are copacetic, unless they aren't |
703 |
|
thisStat = 0 |
704 |
|
|
427 |
– |
!! Fortran's version of a cast: |
428 |
– |
FF_uses_RF = use_RF_c |
429 |
– |
|
705 |
|
!! init_FF is called *after* all of the atom types have been |
706 |
|
!! defined in atype_module using the new_atype subroutine. |
707 |
|
!! |
709 |
|
!! interactions are used by the force field. |
710 |
|
|
711 |
|
FF_uses_DirectionalAtoms = .false. |
437 |
– |
FF_uses_LennardJones = .false. |
438 |
– |
FF_uses_Electrostatics = .false. |
439 |
– |
FF_uses_Charges = .false. |
712 |
|
FF_uses_Dipoles = .false. |
441 |
– |
FF_uses_Sticky = .false. |
442 |
– |
FF_uses_StickyPower = .false. |
713 |
|
FF_uses_GayBerne = .false. |
714 |
|
FF_uses_EAM = .false. |
445 |
– |
FF_uses_Shapes = .false. |
446 |
– |
FF_uses_FLARB = .false. |
715 |
|
|
716 |
|
call getMatchingElementList(atypes, "is_Directional", .true., & |
717 |
|
nMatches, MatchList) |
718 |
|
if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true. |
719 |
|
|
452 |
– |
call getMatchingElementList(atypes, "is_LennardJones", .true., & |
453 |
– |
nMatches, MatchList) |
454 |
– |
if (nMatches .gt. 0) FF_uses_LennardJones = .true. |
455 |
– |
|
456 |
– |
call getMatchingElementList(atypes, "is_Electrostatic", .true., & |
457 |
– |
nMatches, MatchList) |
458 |
– |
if (nMatches .gt. 0) then |
459 |
– |
FF_uses_Electrostatics = .true. |
460 |
– |
endif |
461 |
– |
|
462 |
– |
call getMatchingElementList(atypes, "is_Charge", .true., & |
463 |
– |
nMatches, MatchList) |
464 |
– |
if (nMatches .gt. 0) then |
465 |
– |
FF_uses_Charges = .true. |
466 |
– |
FF_uses_Electrostatics = .true. |
467 |
– |
endif |
468 |
– |
|
720 |
|
call getMatchingElementList(atypes, "is_Dipole", .true., & |
721 |
|
nMatches, MatchList) |
722 |
< |
if (nMatches .gt. 0) then |
472 |
< |
FF_uses_Dipoles = .true. |
473 |
< |
FF_uses_Electrostatics = .true. |
474 |
< |
FF_uses_DirectionalAtoms = .true. |
475 |
< |
endif |
476 |
< |
|
477 |
< |
call getMatchingElementList(atypes, "is_Quadrupole", .true., & |
478 |
< |
nMatches, MatchList) |
479 |
< |
if (nMatches .gt. 0) then |
480 |
< |
FF_uses_Quadrupoles = .true. |
481 |
< |
FF_uses_Electrostatics = .true. |
482 |
< |
FF_uses_DirectionalAtoms = .true. |
483 |
< |
endif |
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 |
722 |
> |
if (nMatches .gt. 0) FF_uses_Dipoles = .true. |
723 |
|
|
724 |
|
call getMatchingElementList(atypes, "is_GayBerne", .true., & |
725 |
|
nMatches, MatchList) |
726 |
< |
if (nMatches .gt. 0) then |
502 |
< |
FF_uses_GayBerne = .true. |
503 |
< |
FF_uses_DirectionalAtoms = .true. |
504 |
< |
endif |
726 |
> |
if (nMatches .gt. 0) FF_uses_GayBerne = .true. |
727 |
|
|
728 |
|
call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList) |
729 |
|
if (nMatches .gt. 0) FF_uses_EAM = .true. |
730 |
|
|
509 |
– |
call getMatchingElementList(atypes, "is_Shape", .true., & |
510 |
– |
nMatches, MatchList) |
511 |
– |
if (nMatches .gt. 0) then |
512 |
– |
FF_uses_Shapes = .true. |
513 |
– |
FF_uses_DirectionalAtoms = .true. |
514 |
– |
endif |
731 |
|
|
516 |
– |
call getMatchingElementList(atypes, "is_FLARB", .true., & |
517 |
– |
nMatches, MatchList) |
518 |
– |
if (nMatches .gt. 0) FF_uses_FLARB = .true. |
519 |
– |
|
520 |
– |
!! Assume sanity (for the sake of argument) |
732 |
|
haveSaneForceField = .true. |
522 |
– |
|
523 |
– |
!! check to make sure the FF_uses_RF setting makes sense |
524 |
– |
|
525 |
– |
if (FF_uses_dipoles) then |
526 |
– |
if (FF_uses_RF) then |
527 |
– |
dielect = getDielect() |
528 |
– |
call initialize_rf(dielect) |
529 |
– |
endif |
530 |
– |
else |
531 |
– |
if (FF_uses_RF) then |
532 |
– |
write(default_error,*) 'Using Reaction Field with no dipoles? Huh?' |
533 |
– |
thisStat = -1 |
534 |
– |
haveSaneForceField = .false. |
535 |
– |
return |
536 |
– |
endif |
537 |
– |
endif |
538 |
– |
|
539 |
– |
!sticky module does not contain check_sticky_FF anymore |
540 |
– |
!if (FF_uses_sticky) then |
541 |
– |
! call check_sticky_FF(my_status) |
542 |
– |
! if (my_status /= 0) then |
543 |
– |
! thisStat = -1 |
544 |
– |
! haveSaneForceField = .false. |
545 |
– |
! return |
546 |
– |
! end if |
547 |
– |
!endif |
733 |
|
|
734 |
|
if (FF_uses_EAM) then |
735 |
|
call init_EAM_FF(my_status) |
741 |
|
end if |
742 |
|
endif |
743 |
|
|
559 |
– |
if (FF_uses_GayBerne) then |
560 |
– |
call check_gb_pair_FF(my_status) |
561 |
– |
if (my_status .ne. 0) then |
562 |
– |
thisStat = -1 |
563 |
– |
haveSaneForceField = .false. |
564 |
– |
return |
565 |
– |
endif |
566 |
– |
endif |
567 |
– |
|
568 |
– |
if (FF_uses_GayBerne .and. FF_uses_LennardJones) then |
569 |
– |
endif |
570 |
– |
|
744 |
|
if (.not. haveNeighborList) then |
745 |
|
!! Create neighbor lists |
746 |
|
call expandNeighborList(nLocal, my_status) |
774 |
|
|
775 |
|
!! Stress Tensor |
776 |
|
real( kind = dp), dimension(9) :: tau |
777 |
< |
real ( kind = dp ) :: pot |
777 |
> |
real ( kind = dp ),dimension(LR_POT_TYPES) :: pot |
778 |
|
logical ( kind = 2) :: do_pot_c, do_stress_c |
779 |
|
logical :: do_pot |
780 |
|
logical :: do_stress |
781 |
|
logical :: in_switching_region |
782 |
|
#ifdef IS_MPI |
783 |
< |
real( kind = DP ) :: pot_local |
783 |
> |
real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local |
784 |
|
integer :: nAtomsInRow |
785 |
|
integer :: nAtomsInCol |
786 |
|
integer :: nprocs |
795 |
|
integer :: nlist |
796 |
|
real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij |
797 |
|
real( kind = DP ) :: sw, dswdr, swderiv, mf |
798 |
+ |
real( kind = DP ) :: rVal |
799 |
|
real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij |
800 |
|
real(kind=dp) :: rfpot, mu_i, virial |
801 |
+ |
real(kind=dp):: rCut |
802 |
|
integer :: me_i, me_j, n_in_i, n_in_j |
803 |
|
logical :: is_dp_i |
804 |
|
integer :: neighborListSize |
806 |
|
integer :: localError |
807 |
|
integer :: propPack_i, propPack_j |
808 |
|
integer :: loopStart, loopEnd, loop |
809 |
< |
integer :: iMap |
810 |
< |
real(kind=dp) :: listSkin = 1.0 |
809 |
> |
integer :: iHash |
810 |
> |
integer :: i1 |
811 |
> |
|
812 |
|
|
813 |
|
!! initialize local variables |
814 |
|
|
872 |
|
! (but only on the first time through): |
873 |
|
if (loop .eq. loopStart) then |
874 |
|
#ifdef IS_MPI |
875 |
< |
call checkNeighborList(nGroupsInRow, q_group_row, listSkin, & |
875 |
> |
call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, & |
876 |
|
update_nlist) |
877 |
|
#else |
878 |
< |
call checkNeighborList(nGroups, q_group, listSkin, & |
878 |
> |
call checkNeighborList(nGroups, q_group, skinThickness, & |
879 |
|
update_nlist) |
880 |
|
#endif |
881 |
|
endif |
926 |
|
endif |
927 |
|
|
928 |
|
#ifdef IS_MPI |
929 |
+ |
me_j = atid_col(j) |
930 |
|
call get_interatomic_vector(q_group_Row(:,i), & |
931 |
|
q_group_Col(:,j), d_grp, rgrpsq) |
932 |
|
#else |
933 |
+ |
me_j = atid(j) |
934 |
|
call get_interatomic_vector(q_group(:,i), & |
935 |
|
q_group(:,j), d_grp, rgrpsq) |
936 |
< |
#endif |
936 |
> |
#endif |
937 |
|
|
938 |
< |
if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then |
938 |
> |
if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then |
939 |
|
if (update_nlist) then |
940 |
|
nlist = nlist + 1 |
941 |
|
|
955 |
|
|
956 |
|
list(nlist) = j |
957 |
|
endif |
958 |
+ |
|
959 |
|
|
960 |
< |
if (loop .eq. PAIR_LOOP) then |
961 |
< |
vij = 0.0d0 |
783 |
< |
fij(1:3) = 0.0d0 |
784 |
< |
endif |
960 |
> |
|
961 |
> |
if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then |
962 |
|
|
963 |
< |
call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, & |
964 |
< |
in_switching_region) |
965 |
< |
|
966 |
< |
n_in_j = groupStartCol(j+1) - groupStartCol(j) |
967 |
< |
|
968 |
< |
do ia = groupStartRow(i), groupStartRow(i+1)-1 |
969 |
< |
|
970 |
< |
atom1 = groupListRow(ia) |
971 |
< |
|
972 |
< |
inner: do jb = groupStartCol(j), groupStartCol(j+1)-1 |
973 |
< |
|
974 |
< |
atom2 = groupListCol(jb) |
975 |
< |
|
976 |
< |
if (skipThisPair(atom1, atom2)) cycle inner |
977 |
< |
|
978 |
< |
if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then |
979 |
< |
d_atm(1:3) = d_grp(1:3) |
980 |
< |
ratmsq = rgrpsq |
981 |
< |
else |
982 |
< |
#ifdef IS_MPI |
983 |
< |
call get_interatomic_vector(q_Row(:,atom1), & |
984 |
< |
q_Col(:,atom2), d_atm, ratmsq) |
985 |
< |
#else |
986 |
< |
call get_interatomic_vector(q(:,atom1), & |
987 |
< |
q(:,atom2), d_atm, ratmsq) |
988 |
< |
#endif |
989 |
< |
endif |
990 |
< |
|
991 |
< |
if (loop .eq. PREPAIR_LOOP) then |
963 |
> |
rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut |
964 |
> |
if (loop .eq. PAIR_LOOP) then |
965 |
> |
vij = 0.0d0 |
966 |
> |
fij(1:3) = 0.0d0 |
967 |
> |
endif |
968 |
> |
|
969 |
> |
call get_switch(rgrpsq, sw, dswdr, rgrp, & |
970 |
> |
group_switch, in_switching_region) |
971 |
> |
|
972 |
> |
n_in_j = groupStartCol(j+1) - groupStartCol(j) |
973 |
> |
|
974 |
> |
do ia = groupStartRow(i), groupStartRow(i+1)-1 |
975 |
> |
|
976 |
> |
atom1 = groupListRow(ia) |
977 |
> |
|
978 |
> |
inner: do jb = groupStartCol(j), groupStartCol(j+1)-1 |
979 |
> |
|
980 |
> |
atom2 = groupListCol(jb) |
981 |
> |
|
982 |
> |
if (skipThisPair(atom1, atom2)) cycle inner |
983 |
> |
|
984 |
> |
if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then |
985 |
> |
d_atm(1:3) = d_grp(1:3) |
986 |
> |
ratmsq = rgrpsq |
987 |
> |
else |
988 |
> |
#ifdef IS_MPI |
989 |
> |
call get_interatomic_vector(q_Row(:,atom1), & |
990 |
> |
q_Col(:,atom2), d_atm, ratmsq) |
991 |
> |
#else |
992 |
> |
call get_interatomic_vector(q(:,atom1), & |
993 |
> |
q(:,atom2), d_atm, ratmsq) |
994 |
> |
#endif |
995 |
> |
endif |
996 |
> |
|
997 |
> |
if (loop .eq. PREPAIR_LOOP) then |
998 |
|
#ifdef IS_MPI |
999 |
< |
call do_prepair(atom1, atom2, ratmsq, d_atm, sw, & |
1000 |
< |
rgrpsq, d_grp, do_pot, do_stress, & |
1001 |
< |
eFrame, A, f, t, pot_local) |
999 |
> |
call do_prepair(atom1, atom2, ratmsq, d_atm, sw, & |
1000 |
> |
rgrpsq, d_grp, rCut, do_pot, do_stress, & |
1001 |
> |
eFrame, A, f, t, pot_local) |
1002 |
|
#else |
1003 |
< |
call do_prepair(atom1, atom2, ratmsq, d_atm, sw, & |
1004 |
< |
rgrpsq, d_grp, do_pot, do_stress, & |
1005 |
< |
eFrame, A, f, t, pot) |
1003 |
> |
call do_prepair(atom1, atom2, ratmsq, d_atm, sw, & |
1004 |
> |
rgrpsq, d_grp, rCut, do_pot, do_stress, & |
1005 |
> |
eFrame, A, f, t, pot) |
1006 |
|
#endif |
1007 |
< |
else |
1007 |
> |
else |
1008 |
|
#ifdef IS_MPI |
1009 |
< |
call do_pair(atom1, atom2, ratmsq, d_atm, sw, & |
1010 |
< |
do_pot, & |
1011 |
< |
eFrame, A, f, t, pot_local, vpair, fpair) |
1009 |
> |
call do_pair(atom1, atom2, ratmsq, d_atm, sw, & |
1010 |
> |
do_pot, eFrame, A, f, t, pot_local, vpair, & |
1011 |
> |
fpair, d_grp, rgrp, rCut) |
1012 |
|
#else |
1013 |
< |
call do_pair(atom1, atom2, ratmsq, d_atm, sw, & |
1014 |
< |
do_pot, & |
1015 |
< |
eFrame, A, f, t, pot, vpair, fpair) |
1013 |
> |
call do_pair(atom1, atom2, ratmsq, d_atm, sw, & |
1014 |
> |
do_pot, eFrame, A, f, t, pot, vpair, fpair, & |
1015 |
> |
d_grp, rgrp, rCut) |
1016 |
|
#endif |
1017 |
+ |
vij = vij + vpair |
1018 |
+ |
fij(1:3) = fij(1:3) + fpair(1:3) |
1019 |
+ |
endif |
1020 |
+ |
enddo inner |
1021 |
+ |
enddo |
1022 |
|
|
1023 |
< |
vij = vij + vpair |
1024 |
< |
fij(1:3) = fij(1:3) + fpair(1:3) |
1025 |
< |
endif |
1026 |
< |
enddo inner |
1027 |
< |
enddo |
1028 |
< |
|
1029 |
< |
if (loop .eq. PAIR_LOOP) then |
1030 |
< |
if (in_switching_region) then |
1031 |
< |
swderiv = vij*dswdr/rgrp |
1032 |
< |
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 |
< |
|
848 |
< |
do ia=groupStartRow(i), groupStartRow(i+1)-1 |
849 |
< |
atom1=groupListRow(ia) |
850 |
< |
mf = mfactRow(atom1) |
1023 |
> |
if (loop .eq. PAIR_LOOP) then |
1024 |
> |
if (in_switching_region) then |
1025 |
> |
swderiv = vij*dswdr/rgrp |
1026 |
> |
fij(1) = fij(1) + swderiv*d_grp(1) |
1027 |
> |
fij(2) = fij(2) + swderiv*d_grp(2) |
1028 |
> |
fij(3) = fij(3) + swderiv*d_grp(3) |
1029 |
> |
|
1030 |
> |
do ia=groupStartRow(i), groupStartRow(i+1)-1 |
1031 |
> |
atom1=groupListRow(ia) |
1032 |
> |
mf = mfactRow(atom1) |
1033 |
|
#ifdef IS_MPI |
1034 |
< |
f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf |
1035 |
< |
f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf |
1036 |
< |
f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf |
1034 |
> |
f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf |
1035 |
> |
f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf |
1036 |
> |
f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf |
1037 |
|
#else |
1038 |
< |
f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf |
1039 |
< |
f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf |
1040 |
< |
f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf |
1038 |
> |
f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf |
1039 |
> |
f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf |
1040 |
> |
f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf |
1041 |
|
#endif |
1042 |
< |
enddo |
1043 |
< |
|
1044 |
< |
do jb=groupStartCol(j), groupStartCol(j+1)-1 |
1045 |
< |
atom2=groupListCol(jb) |
1046 |
< |
mf = mfactCol(atom2) |
1042 |
> |
enddo |
1043 |
> |
|
1044 |
> |
do jb=groupStartCol(j), groupStartCol(j+1)-1 |
1045 |
> |
atom2=groupListCol(jb) |
1046 |
> |
mf = mfactCol(atom2) |
1047 |
|
#ifdef IS_MPI |
1048 |
< |
f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf |
1049 |
< |
f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf |
1050 |
< |
f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf |
1048 |
> |
f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf |
1049 |
> |
f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf |
1050 |
> |
f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf |
1051 |
|
#else |
1052 |
< |
f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf |
1053 |
< |
f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf |
1054 |
< |
f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf |
1052 |
> |
f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf |
1053 |
> |
f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf |
1054 |
> |
f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf |
1055 |
|
#endif |
1056 |
< |
enddo |
1057 |
< |
endif |
1056 |
> |
enddo |
1057 |
> |
endif |
1058 |
|
|
1059 |
< |
if (do_stress) call add_stress_tensor(d_grp, fij) |
1059 |
> |
if (do_stress) call add_stress_tensor(d_grp, fij) |
1060 |
> |
endif |
1061 |
|
endif |
1062 |
< |
end if |
1062 |
> |
endif |
1063 |
|
enddo |
1064 |
+ |
|
1065 |
|
enddo outer |
1066 |
|
|
1067 |
|
if (update_nlist) then |
1121 |
|
|
1122 |
|
if (do_pot) then |
1123 |
|
! scatter/gather pot_row into the members of my column |
1124 |
< |
call scatter(pot_Row, pot_Temp, plan_atom_row) |
1125 |
< |
|
1124 |
> |
do i = 1,LR_POT_TYPES |
1125 |
> |
call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row) |
1126 |
> |
end do |
1127 |
|
! scatter/gather pot_local into all other procs |
1128 |
|
! add resultant to get total pot |
1129 |
|
do i = 1, nlocal |
1130 |
< |
pot_local = pot_local + pot_Temp(i) |
1130 |
> |
pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) & |
1131 |
> |
+ pot_Temp(1:LR_POT_TYPES,i) |
1132 |
|
enddo |
1133 |
|
|
1134 |
|
pot_Temp = 0.0_DP |
1135 |
< |
|
1136 |
< |
call scatter(pot_Col, pot_Temp, plan_atom_col) |
1135 |
> |
do i = 1,LR_POT_TYPES |
1136 |
> |
call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col) |
1137 |
> |
end do |
1138 |
|
do i = 1, nlocal |
1139 |
< |
pot_local = pot_local + pot_Temp(i) |
1139 |
> |
pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)& |
1140 |
> |
+ pot_Temp(1:LR_POT_TYPES,i) |
1141 |
|
enddo |
1142 |
|
|
1143 |
|
endif |
1144 |
|
#endif |
1145 |
|
|
1146 |
< |
if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then |
1146 |
> |
if (SIM_requires_postpair_calc) then |
1147 |
> |
do i = 1, nlocal |
1148 |
> |
|
1149 |
> |
! we loop only over the local atoms, so we don't need row and column |
1150 |
> |
! lookups for the types |
1151 |
> |
|
1152 |
> |
me_i = atid(i) |
1153 |
> |
|
1154 |
> |
! is the atom electrostatic? See if it would have an |
1155 |
> |
! electrostatic interaction with itself |
1156 |
> |
iHash = InteractionHash(me_i,me_i) |
1157 |
|
|
1158 |
< |
if (FF_uses_RF .and. SIM_uses_RF) then |
961 |
< |
|
1158 |
> |
if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then |
1159 |
|
#ifdef IS_MPI |
1160 |
< |
call scatter(rf_Row,rf,plan_atom_row_3d) |
1161 |
< |
call scatter(rf_Col,rf_Temp,plan_atom_col_3d) |
965 |
< |
do i = 1,nlocal |
966 |
< |
rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i) |
967 |
< |
end do |
968 |
< |
#endif |
969 |
< |
|
970 |
< |
do i = 1, nLocal |
971 |
< |
|
972 |
< |
rfpot = 0.0_DP |
973 |
< |
#ifdef IS_MPI |
974 |
< |
me_i = atid_row(i) |
1160 |
> |
call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), & |
1161 |
> |
t, do_pot) |
1162 |
|
#else |
1163 |
< |
me_i = atid(i) |
1163 |
> |
call self_self(i, eFrame, pot(ELECTROSTATIC_POT), & |
1164 |
> |
t, do_pot) |
1165 |
|
#endif |
1166 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
1166 |
> |
endif |
1167 |
> |
|
1168 |
> |
|
1169 |
> |
if (electrostaticSummationMethod.eq.REACTION_FIELD) then |
1170 |
|
|
1171 |
< |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
1172 |
< |
|
1173 |
< |
mu_i = getDipoleMoment(me_i) |
1174 |
< |
|
1175 |
< |
!! The reaction field needs to include a self contribution |
1176 |
< |
!! to the field: |
1177 |
< |
call accumulate_self_rf(i, mu_i, eFrame) |
1178 |
< |
!! Get the reaction field contribution to the |
1179 |
< |
!! potential and torques: |
1180 |
< |
call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot) |
1171 |
> |
! loop over the excludes to accumulate RF stuff we've |
1172 |
> |
! left out of the normal pair loop |
1173 |
> |
|
1174 |
> |
do i1 = 1, nSkipsForAtom(i) |
1175 |
> |
j = skipsForAtom(i, i1) |
1176 |
> |
|
1177 |
> |
! prevent overcounting of the skips |
1178 |
> |
if (i.lt.j) then |
1179 |
> |
call get_interatomic_vector(q(:,i), & |
1180 |
> |
q(:,j), d_atm, ratmsq) |
1181 |
> |
rVal = dsqrt(ratmsq) |
1182 |
> |
call get_switch(ratmsq, sw, dswdr, rVal, group_switch, & |
1183 |
> |
in_switching_region) |
1184 |
|
#ifdef IS_MPI |
1185 |
< |
pot_local = pot_local + rfpot |
1185 |
> |
call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, & |
1186 |
> |
vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot) |
1187 |
|
#else |
1188 |
< |
pot = pot + rfpot |
1189 |
< |
|
1188 |
> |
call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, & |
1189 |
> |
vpair, pot(ELECTROSTATIC_POT), f, t, do_pot) |
1190 |
|
#endif |
1191 |
< |
endif |
1192 |
< |
enddo |
1193 |
< |
endif |
1191 |
> |
endif |
1192 |
> |
enddo |
1193 |
> |
endif |
1194 |
> |
enddo |
1195 |
|
endif |
1196 |
< |
|
1001 |
< |
|
1196 |
> |
|
1197 |
|
#ifdef IS_MPI |
1198 |
< |
|
1198 |
> |
|
1199 |
|
if (do_pot) then |
1200 |
< |
pot = pot + pot_local |
1201 |
< |
!! 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... |
1200 |
> |
call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, & |
1201 |
> |
mpi_comm_world,mpi_err) |
1202 |
|
endif |
1203 |
< |
|
1203 |
> |
|
1204 |
|
if (do_stress) then |
1205 |
|
call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, & |
1206 |
|
mpi_comm_world,mpi_err) |
1207 |
|
call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, & |
1208 |
|
mpi_comm_world,mpi_err) |
1209 |
|
endif |
1210 |
< |
|
1210 |
> |
|
1211 |
|
#else |
1212 |
< |
|
1212 |
> |
|
1213 |
|
if (do_stress) then |
1214 |
|
tau = tau_Temp |
1215 |
|
virial = virial_Temp |
1216 |
|
endif |
1217 |
< |
|
1217 |
> |
|
1218 |
|
#endif |
1219 |
< |
|
1219 |
> |
|
1220 |
|
end subroutine do_force_loop |
1221 |
|
|
1222 |
|
subroutine do_pair(i, j, rijsq, d, sw, do_pot, & |
1223 |
< |
eFrame, A, f, t, pot, vpair, fpair) |
1223 |
> |
eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut) |
1224 |
|
|
1225 |
< |
real( kind = dp ) :: pot, vpair, sw |
1225 |
> |
real( kind = dp ) :: vpair, sw |
1226 |
> |
real( kind = dp ), dimension(LR_POT_TYPES) :: pot |
1227 |
|
real( kind = dp ), dimension(3) :: fpair |
1228 |
|
real( kind = dp ), dimension(nLocal) :: mfact |
1229 |
|
real( kind = dp ), dimension(9,nLocal) :: eFrame |
1234 |
|
logical, intent(inout) :: do_pot |
1235 |
|
integer, intent(in) :: i, j |
1236 |
|
real ( kind = dp ), intent(inout) :: rijsq |
1237 |
< |
real ( kind = dp ) :: r |
1237 |
> |
real ( kind = dp ), intent(inout) :: r_grp |
1238 |
|
real ( kind = dp ), intent(inout) :: d(3) |
1239 |
< |
real ( kind = dp ) :: ebalance |
1239 |
> |
real ( kind = dp ), intent(inout) :: d_grp(3) |
1240 |
> |
real ( kind = dp ), intent(inout) :: rCut |
1241 |
> |
real ( kind = dp ) :: r |
1242 |
|
integer :: me_i, me_j |
1243 |
|
|
1244 |
< |
integer :: iMap |
1244 |
> |
integer :: iHash |
1245 |
|
|
1246 |
|
r = sqrt(rijsq) |
1247 |
|
vpair = 0.0d0 |
1255 |
|
me_j = atid(j) |
1256 |
|
#endif |
1257 |
|
|
1258 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
1259 |
< |
|
1260 |
< |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
1261 |
< |
call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) |
1262 |
< |
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_RF .and. SIM_uses_RF) then |
1072 |
< |
|
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 |
< |
|
1258 |
> |
iHash = InteractionHash(me_i, me_j) |
1259 |
> |
|
1260 |
> |
if ( iand(iHash, LJ_PAIR).ne.0 ) then |
1261 |
> |
call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & |
1262 |
> |
pot(VDW_POT), f, do_pot) |
1263 |
|
endif |
1264 |
< |
|
1265 |
< |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
1264 |
> |
|
1265 |
> |
if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then |
1266 |
> |
call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & |
1267 |
> |
pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot) |
1268 |
> |
endif |
1269 |
> |
|
1270 |
> |
if ( iand(iHash, STICKY_PAIR).ne.0 ) then |
1271 |
|
call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1272 |
< |
pot, A, f, t, do_pot) |
1272 |
> |
pot(HB_POT), A, f, t, do_pot) |
1273 |
|
endif |
1274 |
< |
|
1275 |
< |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
1274 |
> |
|
1275 |
> |
if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then |
1276 |
|
call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1277 |
< |
pot, A, f, t, do_pot) |
1277 |
> |
pot(HB_POT), A, f, t, do_pot) |
1278 |
|
endif |
1279 |
< |
|
1280 |
< |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
1279 |
> |
|
1280 |
> |
if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then |
1281 |
|
call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1282 |
< |
pot, A, f, t, do_pot) |
1282 |
> |
pot(VDW_POT), A, f, t, do_pot) |
1283 |
|
endif |
1284 |
|
|
1285 |
< |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
1286 |
< |
! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1287 |
< |
! pot, A, f, t, do_pot) |
1285 |
> |
if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then |
1286 |
> |
call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & |
1287 |
> |
pot(VDW_POT), A, f, t, do_pot) |
1288 |
|
endif |
1289 |
< |
|
1290 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
1291 |
< |
call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, & |
1292 |
< |
do_pot) |
1289 |
> |
|
1290 |
> |
if ( iand(iHash, EAM_PAIR).ne.0 ) then |
1291 |
> |
call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1292 |
> |
pot(METALLIC_POT), f, do_pot) |
1293 |
|
endif |
1294 |
< |
|
1295 |
< |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
1294 |
> |
|
1295 |
> |
if ( iand(iHash, SHAPE_PAIR).ne.0 ) then |
1296 |
|
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1297 |
< |
pot, A, f, t, do_pot) |
1297 |
> |
pot(VDW_POT), A, f, t, do_pot) |
1298 |
|
endif |
1299 |
< |
|
1300 |
< |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
1299 |
> |
|
1300 |
> |
if ( iand(iHash, SHAPE_LJ).ne.0 ) then |
1301 |
|
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
1302 |
< |
pot, A, f, t, do_pot) |
1302 |
> |
pot(VDW_POT), A, f, t, do_pot) |
1303 |
|
endif |
1304 |
+ |
|
1305 |
+ |
if ( iand(iHash, SC_PAIR).ne.0 ) then |
1306 |
+ |
call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & |
1307 |
+ |
pot(METALLIC_POT), f, do_pot) |
1308 |
+ |
endif |
1309 |
+ |
|
1310 |
|
|
1311 |
+ |
|
1312 |
|
end subroutine do_pair |
1313 |
|
|
1314 |
< |
subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, & |
1314 |
> |
subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, & |
1315 |
|
do_pot, do_stress, eFrame, A, f, t, pot) |
1316 |
|
|
1317 |
< |
real( kind = dp ) :: pot, sw |
1317 |
> |
real( kind = dp ) :: sw |
1318 |
> |
real( kind = dp ), dimension(LR_POT_TYPES) :: pot |
1319 |
|
real( kind = dp ), dimension(9,nLocal) :: eFrame |
1320 |
|
real (kind=dp), dimension(9,nLocal) :: A |
1321 |
|
real (kind=dp), dimension(3,nLocal) :: f |
1323 |
|
|
1324 |
|
logical, intent(inout) :: do_pot, do_stress |
1325 |
|
integer, intent(in) :: i, j |
1326 |
< |
real ( kind = dp ), intent(inout) :: rijsq, rcijsq |
1326 |
> |
real ( kind = dp ), intent(inout) :: rijsq, rcijsq, rCut |
1327 |
|
real ( kind = dp ) :: r, rc |
1328 |
|
real ( kind = dp ), intent(inout) :: d(3), dc(3) |
1329 |
|
|
1330 |
< |
integer :: me_i, me_j, iMap |
1330 |
> |
integer :: me_i, me_j, iHash |
1331 |
|
|
1332 |
+ |
r = sqrt(rijsq) |
1333 |
+ |
|
1334 |
|
#ifdef IS_MPI |
1335 |
|
me_i = atid_row(i) |
1336 |
|
me_j = atid_col(j) |
1339 |
|
me_j = atid(j) |
1340 |
|
#endif |
1341 |
|
|
1342 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
1342 |
> |
iHash = InteractionHash(me_i, me_j) |
1343 |
|
|
1344 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
1345 |
< |
call calc_EAM_prepair_rho(i, j, d, r, rijsq ) |
1344 |
> |
if ( iand(iHash, EAM_PAIR).ne.0 ) then |
1345 |
> |
call calc_EAM_prepair_rho(i, j, d, r, rijsq) |
1346 |
> |
endif |
1347 |
> |
|
1348 |
> |
if ( iand(iHash, SC_PAIR).ne.0 ) then |
1349 |
> |
call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut ) |
1350 |
|
endif |
1351 |
|
|
1352 |
|
end subroutine do_prepair |
1354 |
|
|
1355 |
|
subroutine do_preforce(nlocal,pot) |
1356 |
|
integer :: nlocal |
1357 |
< |
real( kind = dp ) :: pot |
1357 |
> |
real( kind = dp ),dimension(LR_POT_TYPES) :: pot |
1358 |
|
|
1359 |
|
if (FF_uses_EAM .and. SIM_uses_EAM) then |
1360 |
< |
call calc_EAM_preforce_Frho(nlocal,pot) |
1360 |
> |
call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT)) |
1361 |
|
endif |
1362 |
+ |
if (FF_uses_SC .and. SIM_uses_SC) then |
1363 |
+ |
call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT)) |
1364 |
+ |
endif |
1365 |
|
|
1366 |
|
|
1367 |
|
end subroutine do_preforce |
1446 |
|
pot_Col = 0.0_dp |
1447 |
|
pot_Temp = 0.0_dp |
1448 |
|
|
1242 |
– |
rf_Row = 0.0_dp |
1243 |
– |
rf_Col = 0.0_dp |
1244 |
– |
rf_Temp = 0.0_dp |
1245 |
– |
|
1449 |
|
#endif |
1450 |
|
|
1451 |
|
if (FF_uses_EAM .and. SIM_uses_EAM) then |
1452 |
|
call clean_EAM() |
1453 |
|
endif |
1454 |
|
|
1252 |
– |
rf = 0.0_dp |
1455 |
|
tau_Temp = 0.0_dp |
1456 |
|
virial_Temp = 0.0_dp |
1457 |
|
end subroutine zero_work_arrays |
1540 |
|
|
1541 |
|
function FF_UsesDirectionalAtoms() result(doesit) |
1542 |
|
logical :: doesit |
1543 |
< |
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 |
1543 |
> |
doesit = FF_uses_DirectionalAtoms |
1544 |
|
end function FF_UsesDirectionalAtoms |
1545 |
|
|
1546 |
|
function FF_RequiresPrepairCalc() result(doesit) |
1547 |
|
logical :: doesit |
1548 |
< |
doesit = FF_uses_EAM |
1548 |
> |
doesit = FF_uses_EAM .or. FF_uses_SC & |
1549 |
> |
.or. FF_uses_MEAM |
1550 |
|
end function FF_RequiresPrepairCalc |
1551 |
|
|
1351 |
– |
function FF_RequiresPostpairCalc() result(doesit) |
1352 |
– |
logical :: doesit |
1353 |
– |
doesit = FF_uses_RF |
1354 |
– |
end function FF_RequiresPostpairCalc |
1355 |
– |
|
1552 |
|
#ifdef PROFILE |
1553 |
|
function getforcetime() result(totalforcetime) |
1554 |
|
real(kind=dp) :: totalforcetime |