ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE-4/src/UseTheForce/doForces.F90

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2267 by tim, Fri Jul 29 17:34:06 2005 UTC vs.
Revision 2512 by gezelter, Thu Dec 15 21:43:16 2005 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.25 2005-07-29 17:34:06 tim Exp $, $Date: 2005-07-29 17:34:06 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
48	>	!! @version $Id: doForces.F90,v 1.71 2005-12-15 21:43:16 gezelter Exp $, $Date: 2005-12-15 21:43:16 $, $Name: not supported by cvs2svn $, $Revision: 1.71 $
49
50
51		module doForces
#	Line 58 | Line 58 | 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
#	Line 73 | Line 73 | module doForces
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
#	Line 134 | Line 129 | module doForces
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, intent(out) :: status
156		integer :: i
157		integer :: j
158	<	integer :: ihash
157	<	real(kind=dp) :: myRcut
158	>	integer :: iHash
159		!! Test Types
160		logical :: i_is_LJ
161		logical :: i_is_Elect
#	Line 163 | Line 164 | contains
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
#	Line 170 | Line 173 | contains
173		logical :: j_is_GB
174		logical :: j_is_EAM
175		logical :: j_is_Shape
176	<
177	<	status = 0
178	<
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!")
178	<	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)
#	Line 198 | Line 211 | contains
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
#	Line 211 | Line 226 | contains
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)
#	Line 232 | Line 249 | contains
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)
#	Line 241 | Line 262 | contains
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
271
272	<	haveInteractionMap = .true.
273	<	end subroutine createInteractionMap
272	>	haveInteractionHash = .true.
273	>	end subroutine createInteractionHash
274
275	<	! 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.
255	<	subroutine createRcuts(defaultRList,stat)
256	<	real(kind=dp), intent(in), optional :: defaultRList
257	<	integer :: iMap
258	<	integer :: map_i,map_j
259	<	real(kind=dp) :: thisRCut = 0.0_dp
260	<	real(kind=dp) :: actualCutoff = 0.0_dp
261	<	integer, intent(out) :: stat
262	<	integer :: nAtypes
263	<	integer :: myStatus
275	>	subroutine createGtypeCutoffMap()
276
277	<	stat = 0
278	<	if (.not. haveInteractionMap) then
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	<	call createInteractionMap(myStatus)
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 (myStatus .ne. 0) then
295	<	write(default_error, *) 'createInteractionMap failed in doForces!'
272	<	stat = -1
273	<	return
274	<	endif
294	>	if (.not. haveInteractionHash) then
295	>	call createInteractionHash()
296		endif
297	<
298	<
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)
298	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
347	>
348	>	if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
349	>	biggestAtypeCutoff = atypeMaxCutoff(i)
350		endif
300	–
301	–	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
302	–	! thisRCut = getStickyCutOff(map_i,map_j)
303	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
304	–	endif
305	–
306	–	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
307	–	! thisRCut = getStickyPowerCutOff(map_i,map_j)
308	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
309	–	endif
310	–
311	–	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
312	–	! thisRCut = getGayberneCutOff(map_i,map_j)
313	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
314	–	endif
315	–
316	–	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
317	–	!              thisRCut = getGaybrneLJCutOff(map_i,map_j)
318	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
319	–	endif
320	–
321	–	if ( iand(iMap, EAM_PAIR).ne.0 ) then
322	–	!              thisRCut = getEAMCutOff(map_i,map_j)
323	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
324	–	endif
325	–
326	–	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
327	–	!              thisRCut = getShapeCutOff(map_i,map_j)
328	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
329	–	endif
330	–
331	–	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
332	–	!              thisRCut = getShapeLJCutOff(map_i,map_j)
333	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
334	–	endif
335	–	InteractionMap(map_i, map_j)%rList = actualCutoff
336	–	InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
337	–	end do
338	–	end do
339	–	haveRlist = .true.
340	–	end subroutine createRcuts
351
352	<
353	<	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
354	<	!!$  subroutine setRlistDF( this_rlist )
355	<	!!$
356	<	!!$   real(kind=dp) :: this_rlist
357	<	!!$
358	<	!!$    rlist = this_rlist
359	<	!!$    rlistsq = rlist * rlist
360	<	!!$
361	<	!!$    haveRlist = .true.
362	<	!!$
363	<	!!$  end subroutine setRlistDF
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	+	if(.not.allocated(gtypeMaxCutoffRow)) then
381	+	allocate(gtypeMaxCutoffRow(iend))
382	+	else
383	+	deallocate(gtypeMaxCutoffRow)
384	+	allocate(gtypeMaxCutoffRow(iend))
385	+	endif
386
356	–	subroutine setSimVariables()
357	–	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
358	–	SIM_uses_LennardJones = SimUsesLennardJones()
359	–	SIM_uses_Electrostatics = SimUsesElectrostatics()
360	–	SIM_uses_Charges = SimUsesCharges()
361	–	SIM_uses_Dipoles = SimUsesDipoles()
362	–	SIM_uses_Sticky = SimUsesSticky()
363	–	SIM_uses_StickyPower = SimUsesStickyPower()
364	–	SIM_uses_GayBerne = SimUsesGayBerne()
365	–	SIM_uses_EAM = SimUsesEAM()
366	–	SIM_uses_Shapes = SimUsesShapes()
367	–	SIM_uses_FLARB = SimUsesFLARB()
368	–	SIM_uses_RF = SimUsesRF()
369	–	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
370	–	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
371	–	SIM_uses_PBC = SimUsesPBC()
387
388	<	haveSIMvariables = .true.
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	<	return
398	<	end subroutine setSimVariables
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	+	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	+	if (nGroupTypesRow.eq.0) then
440	+	nGroupTypesRow = nGroupTypesRow + 1
441	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
442	+	groupToGtypeRow(i) = nGroupTypesRow
443	+	else
444	+	GtypeFound = .false.
445	+	do g = 1, nGroupTypesRow
446	+	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
447	+	groupToGtypeRow(i) = g
448	+	GtypeFound = .true.
449	+	endif
450	+	enddo
451	+	if (.not.GtypeFound) then
452	+	nGroupTypesRow = nGroupTypesRow + 1
453	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
454	+	groupToGtypeRow(i) = nGroupTypesRow
455	+	endif
456	+	endif
457	+	enddo
458	+
459	+	#ifdef IS_MPI
460	+	do j = jstart, jend
461	+	n_in_j = groupStartCol(j+1) - groupStartCol(j)
462	+	groupMaxCutoffCol(j) = 0.0_dp
463	+	do ja = groupStartCol(j), groupStartCol(j+1)-1
464	+	atom1 = groupListCol(ja)
465	+
466	+	me_j = atid_col(atom1)
467	+
468	+	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
469	+	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
470	+	endif
471	+	enddo
472	+
473	+	if (nGroupTypesCol.eq.0) then
474	+	nGroupTypesCol = nGroupTypesCol + 1
475	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
476	+	groupToGtypeCol(j) = nGroupTypesCol
477	+	else
478	+	GtypeFound = .false.
479	+	do g = 1, nGroupTypesCol
480	+	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
481	+	groupToGtypeCol(j) = g
482	+	GtypeFound = .true.
483	+	endif
484	+	enddo
485	+	if (.not.GtypeFound) then
486	+	nGroupTypesCol = nGroupTypesCol + 1
487	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
488	+	groupToGtypeCol(j) = nGroupTypesCol
489	+	endif
490	+	endif
491	+	enddo
492	+
493	+	#else
494	+	! Set pointers to information we just found
495	+	nGroupTypesCol = nGroupTypesRow
496	+	groupToGtypeCol => groupToGtypeRow
497	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
498	+	groupMaxCutoffCol => groupMaxCutoffRow
499	+	#endif
500	+
501	+	!! allocate the gtypeCutoffMap here.
502	+	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
503	+	!! then we do a double loop over all the group TYPES to find the cutoff
504	+	!! map between groups of two types
505	+	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
506	+
507	+	do i = 1, nGroupTypesRow
508	+	do j = 1, nGroupTypesCol
509	+
510	+	select case(cutoffPolicy)
511	+	case(TRADITIONAL_CUTOFF_POLICY)
512	+	thisRcut = tradRcut
513	+	case(MIX_CUTOFF_POLICY)
514	+	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
515	+	case(MAX_CUTOFF_POLICY)
516	+	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
517	+	case default
518	+	call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
519	+	return
520	+	end select
521	+	gtypeCutoffMap(i,j)%rcut = thisRcut
522	+
523	+	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
524	+
525	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
526	+
527	+	if (.not.haveSkinThickness) then
528	+	skinThickness = 1.0_dp
529	+	endif
530	+
531	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
532	+
533	+	! sanity check
534	+
535	+	if (haveDefaultCutoffs) then
536	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
537	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
538	+	endif
539	+	endif
540	+	enddo
541	+	enddo
542	+
543	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
544	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
545	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
546	+	#ifdef IS_MPI
547	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
548	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
549	+	#endif
550	+	groupMaxCutoffCol => null()
551	+	gtypeMaxCutoffCol => null()
552	+
553	+	haveGtypeCutoffMap = .true.
554	+	end subroutine createGtypeCutoffMap
555	+
556	+	subroutine setCutoffs(defRcut, defRsw)
557	+
558	+	real(kind=dp),intent(in) :: defRcut, defRsw
559	+	character(len = statusMsgSize) :: errMsg
560	+	integer :: localError
561	+
562	+	defaultRcut = defRcut
563	+	defaultRsw = defRsw
564	+
565	+	defaultDoShift = .false.
566	+	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
567	+
568	+	write(errMsg, *) &
569	+	'cutoffRadius and switchingRadius are set to the same', newline &
570	+	// tab, 'value.  OOPSE will use shifted ', newline &
571	+	// tab, 'potentials instead of switching functions.'
572	+
573	+	call handleInfo("setCutoffs", errMsg)
574	+
575	+	defaultDoShift = .true.
576	+
577	+	endif
578	+
579	+	localError = 0
580	+	call setLJDefaultCutoff( defaultRcut, defaultDoShift )
581	+	call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
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	+	haveGtypeCutoffMap = .false.
591	+	end subroutine setCutoffs
592	+
593	+	subroutine cWasLame()
594	+
595	+	VisitCutoffsAfterComputing = .true.
596	+	return
597	+
598	+	end subroutine cWasLame
599	+
600	+	subroutine setCutoffPolicy(cutPolicy)
601	+
602	+	integer, intent(in) :: cutPolicy
603	+
604	+	cutoffPolicy = cutPolicy
605	+	haveCutoffPolicy = .true.
606	+	haveGtypeCutoffMap = .false.
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	+	haveGtypeCutoffMap = .false.
626	+
627	+	end subroutine setSkinThickness
628	+
629	+	subroutine setSimVariables()
630	+	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
631	+	SIM_uses_EAM = SimUsesEAM()
632	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
633	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
634	+	SIM_uses_PBC = SimUsesPBC()
635	+
636	+	haveSIMvariables = .true.
637	+
638	+	return
639	+	end subroutine setSimVariables
640	+
641		subroutine doReadyCheck(error)
642		integer, intent(out) :: error
643
#	Line 382 | Line 645 | contains
645
646		error = 0
647
648	<	if (.not. haveInteractionMap) then
648	>	if (.not. haveInteractionHash) then
649	>	call createInteractionHash()
650	>	endif
651
652	<	myStatus = 0
652	>	if (.not. haveGtypeCutoffMap) then
653	>	call createGtypeCutoffMap()
654	>	endif
655
389	–	call createInteractionMap(myStatus)
656
657	<	if (myStatus .ne. 0) then
658	<	write(default_error, *) 'createInteractionMap failed in doForces!'
393	<	error = -1
394	<	return
395	<	endif
657	>	if (VisitCutoffsAfterComputing) then
658	>	call set_switch(GROUP_SWITCH, largestRcut, largestRcut)
659		endif
660
661	+
662		if (.not. haveSIMvariables) then
663		call setSimVariables()
664		endif
665
666	<	if (.not. haveRlist) then
667	<	write(default_error, *) 'rList has not been set in doForces!'
668	<	error = -1
669	<	return
670	<	endif
666	>	!  if (.not. haveRlist) then
667	>	!     write(default_error, *) 'rList has not been set in doForces!'
668	>	!     error = -1
669	>	!     return
670	>	!  endif
671
672		if (.not. haveNeighborList) then
673		write(default_error, *) 'neighbor list has not been initialized in doForces!'
#	Line 428 | Line 692 | contains
692		end subroutine doReadyCheck
693
694
695	<	subroutine init_FF(use_RF_c, thisStat)
695	>	subroutine init_FF(thisStat)
696
433	–	logical, intent(in) :: use_RF_c
434	–
697		integer, intent(out) :: thisStat
698		integer :: my_status, nMatches
699		integer, pointer :: MatchList(:) => null()
438	–	real(kind=dp) :: rcut, rrf, rt, dielect
700
701		!! assume things are copacetic, unless they aren't
702		thisStat = 0
703
443	–	!! Fortran's version of a cast:
444	–	FF_uses_RF = use_RF_c
445	–
704		!! init_FF is called *after* all of the atom types have been
705		!! defined in atype_module using the new_atype subroutine.
706		!!
#	Line 450 | Line 708 | contains
708		!! interactions are used by the force field.
709
710		FF_uses_DirectionalAtoms = .false.
453	–	FF_uses_LennardJones = .false.
454	–	FF_uses_Electrostatics = .false.
455	–	FF_uses_Charges = .false.
711		FF_uses_Dipoles = .false.
457	–	FF_uses_Sticky = .false.
458	–	FF_uses_StickyPower = .false.
712		FF_uses_GayBerne = .false.
713		FF_uses_EAM = .false.
461	–	FF_uses_Shapes = .false.
462	–	FF_uses_FLARB = .false.
714
715		call getMatchingElementList(atypes, "is_Directional", .true., &
716		nMatches, MatchList)
717		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
718
468	–	call getMatchingElementList(atypes, "is_LennardJones", .true., &
469	–	nMatches, MatchList)
470	–	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
471	–
472	–	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
473	–	nMatches, MatchList)
474	–	if (nMatches .gt. 0) then
475	–	FF_uses_Electrostatics = .true.
476	–	endif
477	–
478	–	call getMatchingElementList(atypes, "is_Charge", .true., &
479	–	nMatches, MatchList)
480	–	if (nMatches .gt. 0) then
481	–	FF_uses_Charges = .true.
482	–	FF_uses_Electrostatics = .true.
483	–	endif
484	–
719		call getMatchingElementList(atypes, "is_Dipole", .true., &
720		nMatches, MatchList)
721	<	if (nMatches .gt. 0) then
488	<	FF_uses_Dipoles = .true.
489	<	FF_uses_Electrostatics = .true.
490	<	FF_uses_DirectionalAtoms = .true.
491	<	endif
492	<
493	<	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
494	<	nMatches, MatchList)
495	<	if (nMatches .gt. 0) then
496	<	FF_uses_Quadrupoles = .true.
497	<	FF_uses_Electrostatics = .true.
498	<	FF_uses_DirectionalAtoms = .true.
499	<	endif
500	<
501	<	call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
502	<	MatchList)
503	<	if (nMatches .gt. 0) then
504	<	FF_uses_Sticky = .true.
505	<	FF_uses_DirectionalAtoms = .true.
506	<	endif
507	<
508	<	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
509	<	MatchList)
510	<	if (nMatches .gt. 0) then
511	<	FF_uses_StickyPower = .true.
512	<	FF_uses_DirectionalAtoms = .true.
513	<	endif
721	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
722
723		call getMatchingElementList(atypes, "is_GayBerne", .true., &
724		nMatches, MatchList)
725	<	if (nMatches .gt. 0) then
518	<	FF_uses_GayBerne = .true.
519	<	FF_uses_DirectionalAtoms = .true.
520	<	endif
725	>	if (nMatches .gt. 0) FF_uses_GayBerne = .true.
726
727		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
728		if (nMatches .gt. 0) FF_uses_EAM = .true.
729
525	–	call getMatchingElementList(atypes, "is_Shape", .true., &
526	–	nMatches, MatchList)
527	–	if (nMatches .gt. 0) then
528	–	FF_uses_Shapes = .true.
529	–	FF_uses_DirectionalAtoms = .true.
530	–	endif
730
532	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
533	–	nMatches, MatchList)
534	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
535	–
536	–	!! Assume sanity (for the sake of argument)
731		haveSaneForceField = .true.
732
539	–	!! check to make sure the FF_uses_RF setting makes sense
540	–
541	–	if (FF_uses_dipoles) then
542	–	if (FF_uses_RF) then
543	–	dielect = getDielect()
544	–	call initialize_rf(dielect)
545	–	endif
546	–	else
547	–	if (FF_uses_RF) then
548	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
549	–	thisStat = -1
550	–	haveSaneForceField = .false.
551	–	return
552	–	endif
553	–	endif
554	–
555	–	!sticky module does not contain check_sticky_FF anymore
556	–	!if (FF_uses_sticky) then
557	–	!   call check_sticky_FF(my_status)
558	–	!   if (my_status /= 0) then
559	–	!      thisStat = -1
560	–	!      haveSaneForceField = .false.
561	–	!      return
562	–	!   end if
563	–	!endif
564	–
733		if (FF_uses_EAM) then
734		call init_EAM_FF(my_status)
735		if (my_status /= 0) then
#	Line 572 | Line 740 | contains
740		end if
741		endif
742
575	–	if (FF_uses_GayBerne) then
576	–	call check_gb_pair_FF(my_status)
577	–	if (my_status .ne. 0) then
578	–	thisStat = -1
579	–	haveSaneForceField = .false.
580	–	return
581	–	endif
582	–	endif
583	–
584	–	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
585	–	endif
586	–
743		if (.not. haveNeighborList) then
744		!! Create neighbor lists
745		call expandNeighborList(nLocal, my_status)
#	Line 617 | Line 773 | contains
773
774		!! Stress Tensor
775		real( kind = dp), dimension(9) :: tau
776	<	real ( kind = dp ) :: pot
776	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
777		logical ( kind = 2) :: do_pot_c, do_stress_c
778		logical :: do_pot
779		logical :: do_stress
780		logical :: in_switching_region
781		#ifdef IS_MPI
782	<	real( kind = DP ) :: pot_local
782	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
783		integer :: nAtomsInRow
784		integer :: nAtomsInCol
785		integer :: nprocs
#	Line 638 | Line 794 | contains
794		integer :: nlist
795		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
796		real( kind = DP ) :: sw, dswdr, swderiv, mf
797	+	real( kind = DP ) :: rVal
798		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
799		real(kind=dp) :: rfpot, mu_i, virial
800	+	real(kind=dp):: rCut
801		integer :: me_i, me_j, n_in_i, n_in_j
802		logical :: is_dp_i
803		integer :: neighborListSize
#	Line 647 | Line 805 | contains
805		integer :: localError
806		integer :: propPack_i, propPack_j
807		integer :: loopStart, loopEnd, loop
808	<	integer :: iMap
809	<	real(kind=dp) :: listSkin = 1.0
808	>	integer :: iHash
809	>	integer :: i1
810	>
811
812		!! initialize local variables
813
#	Line 712 | Line 871 | contains
871		! (but only on the first time through):
872		if (loop .eq. loopStart) then
873		#ifdef IS_MPI
874	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
874	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
875		update_nlist)
876		#else
877	<	call checkNeighborList(nGroups, q_group, listSkin, &
877	>	call checkNeighborList(nGroups, q_group, skinThickness, &
878		update_nlist)
879		#endif
880		endif
#	Line 739 | Line 898 | contains
898		#endif
899		outer: do i = istart, iend
900
742	–	#ifdef IS_MPI
743	–	me_i = atid_row(i)
744	–	#else
745	–	me_i = atid(i)
746	–	#endif
747	–
901		if (update_nlist) point(i) = nlist + 1
902
903		n_in_i = groupStartRow(i+1) - groupStartRow(i)
#	Line 779 | Line 932 | contains
932		me_j = atid(j)
933		call get_interatomic_vector(q_group(:,i), &
934		q_group(:,j), d_grp, rgrpsq)
935	<	#endif
935	>	#endif
936
937	<	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
937	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
938		if (update_nlist) then
939		nlist = nlist + 1
940
#	Line 801 | Line 954 | contains
954
955		list(nlist) = j
956		endif
957	+
958
959	<	if (loop .eq. PAIR_LOOP) then
960	<	vij = 0.0d0
807	<	fij(1:3) = 0.0d0
808	<	endif
959	>
960	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
961
962	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
963	<	in_switching_region)
964	<
965	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
966	<
967	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
968	<
969	<	atom1 = groupListRow(ia)
970	<
971	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
972	<
973	<	atom2 = groupListCol(jb)
974	<
975	<	if (skipThisPair(atom1, atom2)) cycle inner
976	<
977	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
978	<	d_atm(1:3) = d_grp(1:3)
979	<	ratmsq = rgrpsq
980	<	else
981	<	#ifdef IS_MPI
982	<	call get_interatomic_vector(q_Row(:,atom1), &
983	<	q_Col(:,atom2), d_atm, ratmsq)
984	<	#else
985	<	call get_interatomic_vector(q(:,atom1), &
986	<	q(:,atom2), d_atm, ratmsq)
987	<	#endif
988	<	endif
989	<
990	<	if (loop .eq. PREPAIR_LOOP) then
962	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
963	>	if (loop .eq. PAIR_LOOP) then
964	>	vij = 0.0d0
965	>	fij(1:3) = 0.0d0
966	>	endif
967	>
968	>	call get_switch(rgrpsq, sw, dswdr, rgrp, &
969	>	group_switch, in_switching_region)
970	>
971	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
972	>
973	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
974	>
975	>	atom1 = groupListRow(ia)
976	>
977	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
978	>
979	>	atom2 = groupListCol(jb)
980	>
981	>	if (skipThisPair(atom1, atom2))  cycle inner
982	>
983	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
984	>	d_atm(1:3) = d_grp(1:3)
985	>	ratmsq = rgrpsq
986	>	else
987	>	#ifdef IS_MPI
988	>	call get_interatomic_vector(q_Row(:,atom1), &
989	>	q_Col(:,atom2), d_atm, ratmsq)
990	>	#else
991	>	call get_interatomic_vector(q(:,atom1), &
992	>	q(:,atom2), d_atm, ratmsq)
993	>	#endif
994	>	endif
995	>
996	>	if (loop .eq. PREPAIR_LOOP) then
997		#ifdef IS_MPI
998	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
999	<	rgrpsq, d_grp, do_pot, do_stress, &
1000	<	eFrame, A, f, t, pot_local)
998	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
999	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1000	>	eFrame, A, f, t, pot_local)
1001		#else
1002	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1003	<	rgrpsq, d_grp, do_pot, do_stress, &
1004	<	eFrame, A, f, t, pot)
1002	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1003	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1004	>	eFrame, A, f, t, pot)
1005		#endif
1006	<	else
1006	>	else
1007		#ifdef IS_MPI
1008	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1009	<	do_pot, &
1010	<	eFrame, A, f, t, pot_local, vpair, fpair)
1008	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1009	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1010	>	fpair, d_grp, rgrp, rCut)
1011		#else
1012	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1013	<	do_pot,  &
1014	<	eFrame, A, f, t, pot, vpair, fpair)
1012	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1013	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1014	>	d_grp, rgrp, rCut)
1015		#endif
1016	+	vij = vij + vpair
1017	+	fij(1:3) = fij(1:3) + fpair(1:3)
1018	+	endif
1019	+	enddo inner
1020	+	enddo
1021
1022	<	vij = vij + vpair
1023	<	fij(1:3) = fij(1:3) + fpair(1:3)
1024	<	endif
1025	<	enddo inner
1026	<	enddo
1027	<
1028	<	if (loop .eq. PAIR_LOOP) then
1029	<	if (in_switching_region) then
1030	<	swderiv = vij*dswdr/rgrp
1031	<	fij(1) = fij(1) + swderiv*d_grp(1)
869	<	fij(2) = fij(2) + swderiv*d_grp(2)
870	<	fij(3) = fij(3) + swderiv*d_grp(3)
871	<
872	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
873	<	atom1=groupListRow(ia)
874	<	mf = mfactRow(atom1)
1022	>	if (loop .eq. PAIR_LOOP) then
1023	>	if (in_switching_region) then
1024	>	swderiv = vij*dswdr/rgrp
1025	>	fij(1) = fij(1) + swderiv*d_grp(1)
1026	>	fij(2) = fij(2) + swderiv*d_grp(2)
1027	>	fij(3) = fij(3) + swderiv*d_grp(3)
1028	>
1029	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1030	>	atom1=groupListRow(ia)
1031	>	mf = mfactRow(atom1)
1032		#ifdef IS_MPI
1033	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1034	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1035	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1033	>	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1034	>	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1035	>	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1036		#else
1037	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1038	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1039	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1037	>	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1038	>	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1039	>	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1040		#endif
1041	<	enddo
1042	<
1043	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1044	<	atom2=groupListCol(jb)
1045	<	mf = mfactCol(atom2)
1041	>	enddo
1042	>
1043	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1044	>	atom2=groupListCol(jb)
1045	>	mf = mfactCol(atom2)
1046		#ifdef IS_MPI
1047	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1048	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1049	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1047	>	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1048	>	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1049	>	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1050		#else
1051	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1052	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1053	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1051	>	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1052	>	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1053	>	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1054		#endif
1055	<	enddo
1056	<	endif
1055	>	enddo
1056	>	endif
1057
1058	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1058	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1059	>	endif
1060		endif
1061	<	end if
1061	>	endif
1062		enddo
1063	+
1064		enddo outer
1065
1066		if (update_nlist) then
#	Line 961 | Line 1120 | contains
1120
1121		if (do_pot) then
1122		! scatter/gather pot_row into the members of my column
1123	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1124	<
1123	>	do i = 1,LR_POT_TYPES
1124	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1125	>	end do
1126		! scatter/gather pot_local into all other procs
1127		! add resultant to get total pot
1128		do i = 1, nlocal
1129	<	pot_local = pot_local + pot_Temp(i)
1129	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1130	>	+ pot_Temp(1:LR_POT_TYPES,i)
1131		enddo
1132
1133		pot_Temp = 0.0_DP
1134	<
1135	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1134	>	do i = 1,LR_POT_TYPES
1135	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1136	>	end do
1137		do i = 1, nlocal
1138	<	pot_local = pot_local + pot_Temp(i)
1138	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1139	>	+ pot_Temp(1:LR_POT_TYPES,i)
1140		enddo
1141
1142		endif
1143		#endif
1144
1145	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1145	>	if (SIM_requires_postpair_calc) then
1146	>	do i = 1, nlocal
1147	>
1148	>	! we loop only over the local atoms, so we don't need row and column
1149	>	! lookups for the types
1150	>
1151	>	me_i = atid(i)
1152	>
1153	>	! is the atom electrostatic?  See if it would have an
1154	>	! electrostatic interaction with itself
1155	>	iHash = InteractionHash(me_i,me_i)
1156
1157	<	if (FF_uses_RF .and. SIM_uses_RF) then
985	<
1157	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1158		#ifdef IS_MPI
1159	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1160	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
989	<	do i = 1,nlocal
990	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
991	<	end do
992	<	#endif
993	<
994	<	do i = 1, nLocal
995	<
996	<	rfpot = 0.0_DP
997	<	#ifdef IS_MPI
998	<	me_i = atid_row(i)
1159	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1160	>	t, do_pot)
1161		#else
1162	<	me_i = atid(i)
1162	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1163	>	t, do_pot)
1164		#endif
1165	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1165	>	endif
1166	>
1167	>
1168	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1169
1170	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1171	<
1172	<	mu_i = getDipoleMoment(me_i)
1173	<
1174	<	!! The reaction field needs to include a self contribution
1175	<	!! to the field:
1176	<	call accumulate_self_rf(i, mu_i, eFrame)
1177	<	!! Get the reaction field contribution to the
1178	<	!! potential and torques:
1179	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1170	>	! loop over the excludes to accumulate RF stuff we've
1171	>	! left out of the normal pair loop
1172	>
1173	>	do i1 = 1, nSkipsForAtom(i)
1174	>	j = skipsForAtom(i, i1)
1175	>
1176	>	! prevent overcounting of the skips
1177	>	if (i.lt.j) then
1178	>	call get_interatomic_vector(q(:,i), &
1179	>	q(:,j), d_atm, ratmsq)
1180	>	rVal = dsqrt(ratmsq)
1181	>	call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1182	>	in_switching_region)
1183		#ifdef IS_MPI
1184	<	pot_local = pot_local + rfpot
1184	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1185	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1186		#else
1187	<	pot = pot + rfpot
1188	<
1187	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1188	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1189		#endif
1190	<	endif
1191	<	enddo
1192	<	endif
1190	>	endif
1191	>	enddo
1192	>	endif
1193	>	enddo
1194		endif
1195	<
1025	<
1195	>
1196		#ifdef IS_MPI
1197	<
1197	>
1198		if (do_pot) then
1199	<	pot = pot + pot_local
1200	<	!! we assume the c code will do the allreduce to get the total potential
1031	<	!! we could do it right here if we needed to...
1199	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1200	>	mpi_comm_world,mpi_err)
1201		endif
1202	<
1202	>
1203		if (do_stress) then
1204		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1205		mpi_comm_world,mpi_err)
1206		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1207		mpi_comm_world,mpi_err)
1208		endif
1209	<
1209	>
1210		#else
1211	<
1211	>
1212		if (do_stress) then
1213		tau = tau_Temp
1214		virial = virial_Temp
1215		endif
1216	<
1216	>
1217		#endif
1218	<
1218	>
1219		end subroutine do_force_loop
1220
1221		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1222	<	eFrame, A, f, t, pot, vpair, fpair)
1222	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1223
1224	<	real( kind = dp ) :: pot, vpair, sw
1224	>	real( kind = dp ) :: vpair, sw
1225	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1226		real( kind = dp ), dimension(3) :: fpair
1227		real( kind = dp ), dimension(nLocal)   :: mfact
1228		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1063 | Line 1233 | contains
1233		logical, intent(inout) :: do_pot
1234		integer, intent(in) :: i, j
1235		real ( kind = dp ), intent(inout) :: rijsq
1236	<	real ( kind = dp )                :: r
1236	>	real ( kind = dp ), intent(inout) :: r_grp
1237		real ( kind = dp ), intent(inout) :: d(3)
1238	<	real ( kind = dp ) :: ebalance
1238	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1239	>	real ( kind = dp ), intent(inout) :: rCut
1240	>	real ( kind = dp ) :: r
1241		integer :: me_i, me_j
1242
1243	<	integer :: iMap
1243	>	integer :: iHash
1244
1245		r = sqrt(rijsq)
1246		vpair = 0.0d0
#	Line 1082 | Line 1254 | contains
1254		me_j = atid(j)
1255		#endif
1256
1257	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1258	<
1259	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1260	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1261	<	endif
1090	<
1091	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1092	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1093	<	pot, eFrame, f, t, do_pot)
1094	<
1095	<	if (FF_uses_RF .and. SIM_uses_RF) then
1096	<
1097	<	! CHECK ME (RF needs to know about all electrostatic types)
1098	<	call accumulate_rf(i, j, r, eFrame, sw)
1099	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1100	<	endif
1101	<
1257	>	iHash = InteractionHash(me_i, me_j)
1258	>
1259	>	if ( iand(iHash, LJ_PAIR).ne.0 ) then
1260	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1261	>	pot(VDW_POT), f, do_pot)
1262		endif
1263	<
1264	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1263	>
1264	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1265	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1266	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1267	>	endif
1268	>
1269	>	if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1270		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1271	<	pot, A, f, t, do_pot)
1271	>	pot(HB_POT), A, f, t, do_pot)
1272		endif
1273	<
1274	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1273	>
1274	>	if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1275		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1276	<	pot, A, f, t, do_pot)
1276	>	pot(HB_POT), A, f, t, do_pot)
1277		endif
1278	<
1279	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1278	>
1279	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1280		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1281	<	pot, A, f, t, do_pot)
1281	>	pot(VDW_POT), A, f, t, do_pot)
1282		endif
1283
1284	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1285	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1286	<	!           pot, A, f, t, do_pot)
1284	>	if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1285	>	call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1286	>	pot(VDW_POT), A, f, t, do_pot)
1287		endif
1288	<
1289	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1290	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1291	<	do_pot)
1288	>
1289	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1290	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1291	>	pot(METALLIC_POT), f, do_pot)
1292		endif
1293	<
1294	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1293	>
1294	>	if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1295		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1296	<	pot, A, f, t, do_pot)
1296	>	pot(VDW_POT), A, f, t, do_pot)
1297		endif
1298	<
1299	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1298	>
1299	>	if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1300		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1301	<	pot, A, f, t, do_pot)
1301	>	pot(VDW_POT), A, f, t, do_pot)
1302		endif
1303	+
1304	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1305	+	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1306	+	pot(METALLIC_POT), f, do_pot)
1307	+	endif
1308	+
1309
1310	+
1311		end subroutine do_pair
1312
1313	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1313	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1314		do_pot, do_stress, eFrame, A, f, t, pot)
1315
1316	<	real( kind = dp ) :: pot, sw
1316	>	real( kind = dp ) :: sw
1317	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1318		real( kind = dp ), dimension(9,nLocal) :: eFrame
1319		real (kind=dp), dimension(9,nLocal) :: A
1320		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1149 | Line 1322 | contains
1322
1323		logical, intent(inout) :: do_pot, do_stress
1324		integer, intent(in) :: i, j
1325	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1325	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1326		real ( kind = dp )                :: r, rc
1327		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1328
1329	<	integer :: me_i, me_j, iMap
1329	>	integer :: me_i, me_j, iHash
1330
1331	+	r = sqrt(rijsq)
1332	+
1333		#ifdef IS_MPI
1334		me_i = atid_row(i)
1335		me_j = atid_col(j)
#	Line 1163 | Line 1338 | contains
1338		me_j = atid(j)
1339		#endif
1340
1341	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1341	>	iHash = InteractionHash(me_i, me_j)
1342
1343	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1344	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1343	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1344	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1345		endif
1346	+
1347	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1348	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1349	+	endif
1350
1351		end subroutine do_prepair
1352
1353
1354		subroutine do_preforce(nlocal,pot)
1355		integer :: nlocal
1356	<	real( kind = dp ) :: pot
1356	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1357
1358		if (FF_uses_EAM .and. SIM_uses_EAM) then
1359	<	call calc_EAM_preforce_Frho(nlocal,pot)
1359	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1360		endif
1361	+	if (FF_uses_SC .and. SIM_uses_SC) then
1362	+	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1363	+	endif
1364
1365
1366		end subroutine do_preforce
#	Line 1263 | Line 1445 | contains
1445		pot_Col = 0.0_dp
1446		pot_Temp = 0.0_dp
1447
1266	–	rf_Row = 0.0_dp
1267	–	rf_Col = 0.0_dp
1268	–	rf_Temp = 0.0_dp
1269	–
1448		#endif
1449
1450		if (FF_uses_EAM .and. SIM_uses_EAM) then
1451		call clean_EAM()
1452		endif
1453
1276	–	rf = 0.0_dp
1454		tau_Temp = 0.0_dp
1455		virial_Temp = 0.0_dp
1456		end subroutine zero_work_arrays
#	Line 1362 | Line 1539 | contains
1539
1540		function FF_UsesDirectionalAtoms() result(doesit)
1541		logical :: doesit
1542	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1366	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1367	<	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1542	>	doesit = FF_uses_DirectionalAtoms
1543		end function FF_UsesDirectionalAtoms
1544
1545		function FF_RequiresPrepairCalc() result(doesit)
1546		logical :: doesit
1547	<	doesit = FF_uses_EAM
1547	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1548	>	.or. FF_uses_MEAM
1549		end function FF_RequiresPrepairCalc
1550
1375	–	function FF_RequiresPostpairCalc() result(doesit)
1376	–	logical :: doesit
1377	–	doesit = FF_uses_RF
1378	–	end function FF_RequiresPostpairCalc
1379	–
1551		#ifdef PROFILE
1552		function getforcetime() result(totalforcetime)
1553		real(kind=dp) :: totalforcetime

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines