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 2266 by chuckv, Thu Jul 28 22:12:45 2005 UTC vs.
Revision 2503 by gezelter, Thu Dec 8 22:04:40 2005 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.24 2005-07-28 22:12:45 chuckv Exp $, $Date: 2005-07-28 22:12:45 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
48	>	!! @version $Id: doForces.F90,v 1.70 2005-12-08 22:04:40 gezelter Exp $, $Date: 2005-12-08 22:04:40 $, $Name: not supported by cvs2svn $, $Revision: 1.70 $
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
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
#	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
250	–	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,stat)
254	<	real(kind=dp), intent(in), optional :: defaultRList
255	<	integer :: iMap
256	<	integer :: map_i,map_j
257	<	real(kind=dp) :: thisRCut = 0.0_dp
258	<	real(kind=dp) :: actualCutoff = 0.0_dp
259	<	integer, intent(out) :: stat
260	<	integer :: nAtypes
261	<	integer :: myStatus
272	>	haveInteractionHash = .true.
273	>	end subroutine createInteractionHash
274
275	<	stat = 0
264	<	if (.not. haveInteractionMap) then
275	>	subroutine createGtypeCutoffMap()
276
277	<	call createInteractionMap(myStatus)
278	<
279	<	if (myStatus .ne. 0) then
280	<	write(default_error, *) 'createInteractionMap failed in doForces!'
281	<	stat = -1
282	<	return
283	<	endif
284	<	endif
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)
296	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
347	>
348	>	if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
349	>	biggestAtypeCutoff = atypeMaxCutoff(i)
350		endif
298	–
299	–	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
300	–	!             thisRCut = getStickyCutOff(map_i,map_j)
301	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
302	–	endif
303	–
304	–	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
305	–	!              thisRCut = getStickyPowerCutOff(map_i,map_j)
306	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
307	–	endif
308	–
309	–	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
310	–	!              thisRCut = getGayberneCutOff(map_i,map_j)
311	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
312	–	endif
313	–
314	–	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
315	–	!              thisRCut = getGaybrneLJCutOff(map_i,map_j)
316	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
317	–	endif
318	–
319	–	if ( iand(iMap, EAM_PAIR).ne.0 ) then
320	–	!              thisRCut = getEAMCutOff(map_i,map_j)
321	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
322	–	endif
323	–
324	–	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
325	–	!              thisRCut = getShapeCutOff(map_i,map_j)
326	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
327	–	endif
328	–
329	–	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
330	–	!              thisRCut = getShapeLJCutOff(map_i,map_j)
331	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
332	–	endif
333	–	InteractionMap(map_i, map_j)%rList = actualCutoff
334	–	InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
335	–	end do
336	–	end do
337	–	haveRlist = .true.
338	–	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
354	–	subroutine setSimVariables()
355	–	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
356	–	SIM_uses_LennardJones = SimUsesLennardJones()
357	–	SIM_uses_Electrostatics = SimUsesElectrostatics()
358	–	SIM_uses_Charges = SimUsesCharges()
359	–	SIM_uses_Dipoles = SimUsesDipoles()
360	–	SIM_uses_Sticky = SimUsesSticky()
361	–	SIM_uses_StickyPower = SimUsesStickyPower()
362	–	SIM_uses_GayBerne = SimUsesGayBerne()
363	–	SIM_uses_EAM = SimUsesEAM()
364	–	SIM_uses_Shapes = SimUsesShapes()
365	–	SIM_uses_FLARB = SimUsesFLARB()
366	–	SIM_uses_RF = SimUsesRF()
367	–	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
368	–	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
369	–	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	+
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	+	write(*,*) 'have cutoffPolicy in F = ', cutPolicy
606	+
607	+	call createGtypeCutoffMap()
608	+
609	+	end subroutine setCutoffPolicy
610	+
611	+	subroutine setElectrostaticMethod( thisESM )
612	+
613	+	integer, intent(in) :: thisESM
614	+
615	+	electrostaticSummationMethod = thisESM
616	+	haveElectrostaticSummationMethod = .true.
617	+
618	+	end subroutine setElectrostaticMethod
619	+
620	+	subroutine setSkinThickness( thisSkin )
621	+
622	+	real(kind=dp), intent(in) :: thisSkin
623	+
624	+	skinThickness = thisSkin
625	+	haveSkinThickness = .true.
626	+
627	+	call createGtypeCutoffMap()
628	+
629	+	end subroutine setSkinThickness
630	+
631	+	subroutine setSimVariables()
632	+	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
633	+	SIM_uses_EAM = SimUsesEAM()
634	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
635	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
636	+	SIM_uses_PBC = SimUsesPBC()
637	+
638	+	haveSIMvariables = .true.
639	+
640	+	return
641	+	end subroutine setSimVariables
642	+
643		subroutine doReadyCheck(error)
644		integer, intent(out) :: error
645
#	Line 380 | Line 647 | contains
647
648		error = 0
649
650	<	if (.not. haveInteractionMap) then
650	>	if (.not. haveInteractionHash) then
651	>	call createInteractionHash()
652	>	endif
653
654	<	myStatus = 0
654	>	if (.not. haveGtypeCutoffMap) then
655	>	call createGtypeCutoffMap()
656	>	endif
657
387	–	call createInteractionMap(myStatus)
658
659	<	if (myStatus .ne. 0) then
660	<	write(default_error, *) 'createInteractionMap failed in doForces!'
391	<	error = -1
392	<	return
393	<	endif
659	>	if (VisitCutoffsAfterComputing) then
660	>	call set_switch(GROUP_SWITCH, largestRcut, largestRcut)
661		endif
662
663	+
664		if (.not. haveSIMvariables) then
665		call setSimVariables()
666		endif
667
668	<	if (.not. haveRlist) then
669	<	write(default_error, *) 'rList has not been set in doForces!'
670	<	error = -1
671	<	return
672	<	endif
668	>	!  if (.not. haveRlist) then
669	>	!     write(default_error, *) 'rList has not been set in doForces!'
670	>	!     error = -1
671	>	!     return
672	>	!  endif
673
674		if (.not. haveNeighborList) then
675		write(default_error, *) 'neighbor list has not been initialized in doForces!'
#	Line 426 | Line 694 | contains
694		end subroutine doReadyCheck
695
696
697	<	subroutine init_FF(use_RF_c, thisStat)
697	>	subroutine init_FF(thisStat)
698
431	–	logical, intent(in) :: use_RF_c
432	–
699		integer, intent(out) :: thisStat
700		integer :: my_status, nMatches
701		integer, pointer :: MatchList(:) => null()
436	–	real(kind=dp) :: rcut, rrf, rt, dielect
702
703		!! assume things are copacetic, unless they aren't
704		thisStat = 0
705
441	–	!! Fortran's version of a cast:
442	–	FF_uses_RF = use_RF_c
443	–
706		!! init_FF is called *after* all of the atom types have been
707		!! defined in atype_module using the new_atype subroutine.
708		!!
#	Line 448 | Line 710 | contains
710		!! interactions are used by the force field.
711
712		FF_uses_DirectionalAtoms = .false.
451	–	FF_uses_LennardJones = .false.
452	–	FF_uses_Electrostatics = .false.
453	–	FF_uses_Charges = .false.
713		FF_uses_Dipoles = .false.
455	–	FF_uses_Sticky = .false.
456	–	FF_uses_StickyPower = .false.
714		FF_uses_GayBerne = .false.
715		FF_uses_EAM = .false.
459	–	FF_uses_Shapes = .false.
460	–	FF_uses_FLARB = .false.
716
717		call getMatchingElementList(atypes, "is_Directional", .true., &
718		nMatches, MatchList)
719		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
720
466	–	call getMatchingElementList(atypes, "is_LennardJones", .true., &
467	–	nMatches, MatchList)
468	–	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
469	–
470	–	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
471	–	nMatches, MatchList)
472	–	if (nMatches .gt. 0) then
473	–	FF_uses_Electrostatics = .true.
474	–	endif
475	–
476	–	call getMatchingElementList(atypes, "is_Charge", .true., &
477	–	nMatches, MatchList)
478	–	if (nMatches .gt. 0) then
479	–	FF_uses_Charges = .true.
480	–	FF_uses_Electrostatics = .true.
481	–	endif
482	–
721		call getMatchingElementList(atypes, "is_Dipole", .true., &
722		nMatches, MatchList)
723	<	if (nMatches .gt. 0) then
486	<	FF_uses_Dipoles = .true.
487	<	FF_uses_Electrostatics = .true.
488	<	FF_uses_DirectionalAtoms = .true.
489	<	endif
490	<
491	<	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
492	<	nMatches, MatchList)
493	<	if (nMatches .gt. 0) then
494	<	FF_uses_Quadrupoles = .true.
495	<	FF_uses_Electrostatics = .true.
496	<	FF_uses_DirectionalAtoms = .true.
497	<	endif
498	<
499	<	call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
500	<	MatchList)
501	<	if (nMatches .gt. 0) then
502	<	FF_uses_Sticky = .true.
503	<	FF_uses_DirectionalAtoms = .true.
504	<	endif
505	<
506	<	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
507	<	MatchList)
508	<	if (nMatches .gt. 0) then
509	<	FF_uses_StickyPower = .true.
510	<	FF_uses_DirectionalAtoms = .true.
511	<	endif
723	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
724
725		call getMatchingElementList(atypes, "is_GayBerne", .true., &
726		nMatches, MatchList)
727	<	if (nMatches .gt. 0) then
516	<	FF_uses_GayBerne = .true.
517	<	FF_uses_DirectionalAtoms = .true.
518	<	endif
727	>	if (nMatches .gt. 0) FF_uses_GayBerne = .true.
728
729		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
730		if (nMatches .gt. 0) FF_uses_EAM = .true.
731
523	–	call getMatchingElementList(atypes, "is_Shape", .true., &
524	–	nMatches, MatchList)
525	–	if (nMatches .gt. 0) then
526	–	FF_uses_Shapes = .true.
527	–	FF_uses_DirectionalAtoms = .true.
528	–	endif
732
530	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
531	–	nMatches, MatchList)
532	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
533	–
534	–	!! Assume sanity (for the sake of argument)
733		haveSaneForceField = .true.
734
537	–	!! check to make sure the FF_uses_RF setting makes sense
538	–
539	–	if (FF_uses_dipoles) then
540	–	if (FF_uses_RF) then
541	–	dielect = getDielect()
542	–	call initialize_rf(dielect)
543	–	endif
544	–	else
545	–	if (FF_uses_RF) then
546	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
547	–	thisStat = -1
548	–	haveSaneForceField = .false.
549	–	return
550	–	endif
551	–	endif
552	–
553	–	!sticky module does not contain check_sticky_FF anymore
554	–	!if (FF_uses_sticky) then
555	–	!   call check_sticky_FF(my_status)
556	–	!   if (my_status /= 0) then
557	–	!      thisStat = -1
558	–	!      haveSaneForceField = .false.
559	–	!      return
560	–	!   end if
561	–	!endif
562	–
735		if (FF_uses_EAM) then
736		call init_EAM_FF(my_status)
737		if (my_status /= 0) then
#	Line 570 | Line 742 | contains
742		end if
743		endif
744
573	–	if (FF_uses_GayBerne) then
574	–	call check_gb_pair_FF(my_status)
575	–	if (my_status .ne. 0) then
576	–	thisStat = -1
577	–	haveSaneForceField = .false.
578	–	return
579	–	endif
580	–	endif
581	–
582	–	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
583	–	endif
584	–
745		if (.not. haveNeighborList) then
746		!! Create neighbor lists
747		call expandNeighborList(nLocal, my_status)
#	Line 615 | Line 775 | contains
775
776		!! Stress Tensor
777		real( kind = dp), dimension(9) :: tau
778	<	real ( kind = dp ) :: pot
778	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
779		logical ( kind = 2) :: do_pot_c, do_stress_c
780		logical :: do_pot
781		logical :: do_stress
782		logical :: in_switching_region
783		#ifdef IS_MPI
784	<	real( kind = DP ) :: pot_local
784	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
785		integer :: nAtomsInRow
786		integer :: nAtomsInCol
787		integer :: nprocs
#	Line 636 | Line 796 | contains
796		integer :: nlist
797		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
798		real( kind = DP ) :: sw, dswdr, swderiv, mf
799	+	real( kind = DP ) :: rVal
800		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
801		real(kind=dp) :: rfpot, mu_i, virial
802	+	real(kind=dp):: rCut
803		integer :: me_i, me_j, n_in_i, n_in_j
804		logical :: is_dp_i
805		integer :: neighborListSize
#	Line 645 | Line 807 | contains
807		integer :: localError
808		integer :: propPack_i, propPack_j
809		integer :: loopStart, loopEnd, loop
810	<	integer :: iMap
811	<	real(kind=dp) :: listSkin = 1.0
810	>	integer :: iHash
811	>	integer :: i1
812	>
813
814		!! initialize local variables
815
#	Line 710 | Line 873 | contains
873		! (but only on the first time through):
874		if (loop .eq. loopStart) then
875		#ifdef IS_MPI
876	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
876	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
877		update_nlist)
878		#else
879	<	call checkNeighborList(nGroups, q_group, listSkin, &
879	>	call checkNeighborList(nGroups, q_group, skinThickness, &
880		update_nlist)
881		#endif
882		endif
#	Line 743 | Line 906 | contains
906
907		if (update_nlist) then
908		#ifdef IS_MPI
746	–	me_i = atid_row(i)
909		jstart = 1
910		jend = nGroupsInCol
911		#else
750	–	me_i = atid(i)
912		jstart = i+1
913		jend = nGroups
914		#endif
#	Line 773 | Line 934 | contains
934		me_j = atid(j)
935		call get_interatomic_vector(q_group(:,i), &
936		q_group(:,j), d_grp, rgrpsq)
937	<	#endif
937	>	#endif
938
939	<	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
939	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
940		if (update_nlist) then
941		nlist = nlist + 1
942
#	Line 795 | Line 956 | contains
956
957		list(nlist) = j
958		endif
959	+
960
961	<	if (loop .eq. PAIR_LOOP) then
962	<	vij = 0.0d0
801	<	fij(1:3) = 0.0d0
802	<	endif
961	>
962	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
963
964	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
965	<	in_switching_region)
966	<
967	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
968	<
969	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
970	<
971	<	atom1 = groupListRow(ia)
972	<
973	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
974	<
975	<	atom2 = groupListCol(jb)
976	<
977	<	if (skipThisPair(atom1, atom2)) cycle inner
978	<
979	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
980	<	d_atm(1:3) = d_grp(1:3)
981	<	ratmsq = rgrpsq
982	<	else
983	<	#ifdef IS_MPI
984	<	call get_interatomic_vector(q_Row(:,atom1), &
985	<	q_Col(:,atom2), d_atm, ratmsq)
986	<	#else
987	<	call get_interatomic_vector(q(:,atom1), &
988	<	q(:,atom2), d_atm, ratmsq)
989	<	#endif
990	<	endif
991	<
992	<	if (loop .eq. PREPAIR_LOOP) then
964	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
965	>	if (loop .eq. PAIR_LOOP) then
966	>	vij = 0.0d0
967	>	fij(1:3) = 0.0d0
968	>	endif
969	>
970	>	call get_switch(rgrpsq, sw, dswdr, rgrp, &
971	>	group_switch, in_switching_region)
972	>
973	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
974	>
975	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
976	>
977	>	atom1 = groupListRow(ia)
978	>
979	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
980	>
981	>	atom2 = groupListCol(jb)
982	>
983	>	if (skipThisPair(atom1, atom2))  cycle inner
984	>
985	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
986	>	d_atm(1:3) = d_grp(1:3)
987	>	ratmsq = rgrpsq
988	>	else
989	>	#ifdef IS_MPI
990	>	call get_interatomic_vector(q_Row(:,atom1), &
991	>	q_Col(:,atom2), d_atm, ratmsq)
992	>	#else
993	>	call get_interatomic_vector(q(:,atom1), &
994	>	q(:,atom2), d_atm, ratmsq)
995	>	#endif
996	>	endif
997	>
998	>	if (loop .eq. PREPAIR_LOOP) then
999		#ifdef IS_MPI
1000	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1001	<	rgrpsq, d_grp, do_pot, do_stress, &
1002	<	eFrame, A, f, t, pot_local)
1000	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1001	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1002	>	eFrame, A, f, t, pot_local)
1003		#else
1004	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1005	<	rgrpsq, d_grp, do_pot, do_stress, &
1006	<	eFrame, A, f, t, pot)
1004	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1005	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1006	>	eFrame, A, f, t, pot)
1007		#endif
1008	<	else
1008	>	else
1009		#ifdef IS_MPI
1010	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1011	<	do_pot, &
1012	<	eFrame, A, f, t, pot_local, vpair, fpair)
1010	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1011	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1012	>	fpair, d_grp, rgrp, rCut)
1013		#else
1014	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1015	<	do_pot,  &
1016	<	eFrame, A, f, t, pot, vpair, fpair)
1014	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1015	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1016	>	d_grp, rgrp, rCut)
1017		#endif
1018	+	vij = vij + vpair
1019	+	fij(1:3) = fij(1:3) + fpair(1:3)
1020	+	endif
1021	+	enddo inner
1022	+	enddo
1023
1024	<	vij = vij + vpair
1025	<	fij(1:3) = fij(1:3) + fpair(1:3)
1026	<	endif
1027	<	enddo inner
1028	<	enddo
1029	<
1030	<	if (loop .eq. PAIR_LOOP) then
1031	<	if (in_switching_region) then
1032	<	swderiv = vij*dswdr/rgrp
1033	<	fij(1) = fij(1) + swderiv*d_grp(1)
863	<	fij(2) = fij(2) + swderiv*d_grp(2)
864	<	fij(3) = fij(3) + swderiv*d_grp(3)
865	<
866	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
867	<	atom1=groupListRow(ia)
868	<	mf = mfactRow(atom1)
1024	>	if (loop .eq. PAIR_LOOP) then
1025	>	if (in_switching_region) then
1026	>	swderiv = vij*dswdr/rgrp
1027	>	fij(1) = fij(1) + swderiv*d_grp(1)
1028	>	fij(2) = fij(2) + swderiv*d_grp(2)
1029	>	fij(3) = fij(3) + swderiv*d_grp(3)
1030	>
1031	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1032	>	atom1=groupListRow(ia)
1033	>	mf = mfactRow(atom1)
1034		#ifdef IS_MPI
1035	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1036	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1037	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1035	>	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1036	>	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1037	>	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1038		#else
1039	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1040	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1041	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1039	>	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1040	>	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1041	>	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1042		#endif
1043	<	enddo
1044	<
1045	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1046	<	atom2=groupListCol(jb)
1047	<	mf = mfactCol(atom2)
1043	>	enddo
1044	>
1045	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1046	>	atom2=groupListCol(jb)
1047	>	mf = mfactCol(atom2)
1048		#ifdef IS_MPI
1049	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1050	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1051	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1049	>	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1050	>	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1051	>	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1052		#else
1053	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1054	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1055	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1053	>	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1054	>	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1055	>	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1056		#endif
1057	<	enddo
1058	<	endif
1057	>	enddo
1058	>	endif
1059
1060	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1060	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1061	>	endif
1062		endif
1063	<	end if
1063	>	endif
1064		enddo
1065	+
1066		enddo outer
1067
1068		if (update_nlist) then
#	Line 955 | Line 1122 | contains
1122
1123		if (do_pot) then
1124		! scatter/gather pot_row into the members of my column
1125	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1126	<
1125	>	do i = 1,LR_POT_TYPES
1126	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1127	>	end do
1128		! scatter/gather pot_local into all other procs
1129		! add resultant to get total pot
1130		do i = 1, nlocal
1131	<	pot_local = pot_local + pot_Temp(i)
1131	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1132	>	+ pot_Temp(1:LR_POT_TYPES,i)
1133		enddo
1134
1135		pot_Temp = 0.0_DP
1136	<
1137	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1136	>	do i = 1,LR_POT_TYPES
1137	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1138	>	end do
1139		do i = 1, nlocal
1140	<	pot_local = pot_local + pot_Temp(i)
1140	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1141	>	+ pot_Temp(1:LR_POT_TYPES,i)
1142		enddo
1143
1144		endif
1145		#endif
1146
1147	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1147	>	if (SIM_requires_postpair_calc) then
1148	>	do i = 1, nlocal
1149	>
1150	>	! we loop only over the local atoms, so we don't need row and column
1151	>	! lookups for the types
1152	>
1153	>	me_i = atid(i)
1154	>
1155	>	! is the atom electrostatic?  See if it would have an
1156	>	! electrostatic interaction with itself
1157	>	iHash = InteractionHash(me_i,me_i)
1158
1159	<	if (FF_uses_RF .and. SIM_uses_RF) then
979	<
1159	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1160		#ifdef IS_MPI
1161	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1162	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
983	<	do i = 1,nlocal
984	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
985	<	end do
986	<	#endif
987	<
988	<	do i = 1, nLocal
989	<
990	<	rfpot = 0.0_DP
991	<	#ifdef IS_MPI
992	<	me_i = atid_row(i)
1161	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1162	>	t, do_pot)
1163		#else
1164	<	me_i = atid(i)
1164	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1165	>	t, do_pot)
1166		#endif
1167	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1167	>	endif
1168	>
1169	>
1170	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1171
1172	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1173	<
1174	<	mu_i = getDipoleMoment(me_i)
1175	<
1176	<	!! The reaction field needs to include a self contribution
1177	<	!! to the field:
1178	<	call accumulate_self_rf(i, mu_i, eFrame)
1179	<	!! Get the reaction field contribution to the
1180	<	!! potential and torques:
1181	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1172	>	! loop over the excludes to accumulate RF stuff we've
1173	>	! left out of the normal pair loop
1174	>
1175	>	do i1 = 1, nSkipsForAtom(i)
1176	>	j = skipsForAtom(i, i1)
1177	>
1178	>	! prevent overcounting of the skips
1179	>	if (i.lt.j) then
1180	>	call get_interatomic_vector(q(:,i), &
1181	>	q(:,j), d_atm, ratmsq)
1182	>	rVal = dsqrt(ratmsq)
1183	>	call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1184	>	in_switching_region)
1185		#ifdef IS_MPI
1186	<	pot_local = pot_local + rfpot
1186	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1187	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1188		#else
1189	<	pot = pot + rfpot
1190	<
1189	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1190	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1191		#endif
1192	<	endif
1193	<	enddo
1194	<	endif
1192	>	endif
1193	>	enddo
1194	>	endif
1195	>	enddo
1196		endif
1197	<
1019	<
1197	>
1198		#ifdef IS_MPI
1199	<
1199	>
1200		if (do_pot) then
1201	<	pot = pot + pot_local
1202	<	!! we assume the c code will do the allreduce to get the total potential
1025	<	!! we could do it right here if we needed to...
1201	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1202	>	mpi_comm_world,mpi_err)
1203		endif
1204	<
1204	>
1205		if (do_stress) then
1206		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1207		mpi_comm_world,mpi_err)
1208		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1209		mpi_comm_world,mpi_err)
1210		endif
1211	<
1211	>
1212		#else
1213	<
1213	>
1214		if (do_stress) then
1215		tau = tau_Temp
1216		virial = virial_Temp
1217		endif
1218	<
1218	>
1219		#endif
1220	<
1220	>
1221		end subroutine do_force_loop
1222
1223		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1224	<	eFrame, A, f, t, pot, vpair, fpair)
1224	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1225
1226	<	real( kind = dp ) :: pot, vpair, sw
1226	>	real( kind = dp ) :: vpair, sw
1227	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1228		real( kind = dp ), dimension(3) :: fpair
1229		real( kind = dp ), dimension(nLocal)   :: mfact
1230		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1057 | Line 1235 | contains
1235		logical, intent(inout) :: do_pot
1236		integer, intent(in) :: i, j
1237		real ( kind = dp ), intent(inout) :: rijsq
1238	<	real ( kind = dp )                :: r
1238	>	real ( kind = dp ), intent(inout) :: r_grp
1239		real ( kind = dp ), intent(inout) :: d(3)
1240	<	real ( kind = dp ) :: ebalance
1240	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1241	>	real ( kind = dp ), intent(inout) :: rCut
1242	>	real ( kind = dp ) :: r
1243		integer :: me_i, me_j
1244
1245	<	integer :: iMap
1245	>	integer :: iHash
1246
1247		r = sqrt(rijsq)
1248		vpair = 0.0d0
#	Line 1076 | Line 1256 | contains
1256		me_j = atid(j)
1257		#endif
1258
1259	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1260	<
1261	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1262	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1263	<	endif
1084	<
1085	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1086	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1087	<	pot, eFrame, f, t, do_pot)
1088	<
1089	<	if (FF_uses_RF .and. SIM_uses_RF) then
1090	<
1091	<	! CHECK ME (RF needs to know about all electrostatic types)
1092	<	call accumulate_rf(i, j, r, eFrame, sw)
1093	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1094	<	endif
1095	<
1259	>	iHash = InteractionHash(me_i, me_j)
1260	>
1261	>	if ( iand(iHash, LJ_PAIR).ne.0 ) then
1262	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1263	>	pot(VDW_POT), f, do_pot)
1264		endif
1265	<
1266	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1265	>
1266	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1267	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1268	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1269	>	endif
1270	>
1271	>	if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1272		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1273	<	pot, A, f, t, do_pot)
1273	>	pot(HB_POT), A, f, t, do_pot)
1274		endif
1275	<
1276	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1275	>
1276	>	if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1277		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1278	<	pot, A, f, t, do_pot)
1278	>	pot(HB_POT), A, f, t, do_pot)
1279		endif
1280	<
1281	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1280	>
1281	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1282		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1283	<	pot, A, f, t, do_pot)
1283	>	pot(VDW_POT), A, f, t, do_pot)
1284		endif
1285
1286	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1287	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1288	<	!           pot, A, f, t, do_pot)
1286	>	if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1287	>	call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1288	>	pot(VDW_POT), A, f, t, do_pot)
1289		endif
1290	<
1291	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1292	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1293	<	do_pot)
1290	>
1291	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1292	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1293	>	pot(METALLIC_POT), f, do_pot)
1294		endif
1295	<
1296	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1295	>
1296	>	if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1297		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1298	<	pot, A, f, t, do_pot)
1298	>	pot(VDW_POT), A, f, t, do_pot)
1299		endif
1300	<
1301	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1300	>
1301	>	if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1302		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1303	<	pot, A, f, t, do_pot)
1303	>	pot(VDW_POT), A, f, t, do_pot)
1304		endif
1305	+
1306	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1307	+	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1308	+	pot(METALLIC_POT), f, do_pot)
1309	+	endif
1310	+
1311
1312	+
1313		end subroutine do_pair
1314
1315	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1315	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1316		do_pot, do_stress, eFrame, A, f, t, pot)
1317
1318	<	real( kind = dp ) :: pot, sw
1318	>	real( kind = dp ) :: sw
1319	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1320		real( kind = dp ), dimension(9,nLocal) :: eFrame
1321		real (kind=dp), dimension(9,nLocal) :: A
1322		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1143 | Line 1324 | contains
1324
1325		logical, intent(inout) :: do_pot, do_stress
1326		integer, intent(in) :: i, j
1327	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1327	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1328		real ( kind = dp )                :: r, rc
1329		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1330
1331	<	integer :: me_i, me_j, iMap
1331	>	integer :: me_i, me_j, iHash
1332
1333	+	r = sqrt(rijsq)
1334	+
1335		#ifdef IS_MPI
1336		me_i = atid_row(i)
1337		me_j = atid_col(j)
#	Line 1157 | Line 1340 | contains
1340		me_j = atid(j)
1341		#endif
1342
1343	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1343	>	iHash = InteractionHash(me_i, me_j)
1344
1345	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1346	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1345	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1346	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1347		endif
1348	+
1349	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1350	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1351	+	endif
1352
1353		end subroutine do_prepair
1354
1355
1356		subroutine do_preforce(nlocal,pot)
1357		integer :: nlocal
1358	<	real( kind = dp ) :: pot
1358	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1359
1360		if (FF_uses_EAM .and. SIM_uses_EAM) then
1361	<	call calc_EAM_preforce_Frho(nlocal,pot)
1361	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1362		endif
1363	+	if (FF_uses_SC .and. SIM_uses_SC) then
1364	+	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1365	+	endif
1366
1367
1368		end subroutine do_preforce
#	Line 1257 | Line 1447 | contains
1447		pot_Col = 0.0_dp
1448		pot_Temp = 0.0_dp
1449
1260	–	rf_Row = 0.0_dp
1261	–	rf_Col = 0.0_dp
1262	–	rf_Temp = 0.0_dp
1263	–
1450		#endif
1451
1452		if (FF_uses_EAM .and. SIM_uses_EAM) then
1453		call clean_EAM()
1454		endif
1455
1270	–	rf = 0.0_dp
1456		tau_Temp = 0.0_dp
1457		virial_Temp = 0.0_dp
1458		end subroutine zero_work_arrays
#	Line 1356 | Line 1541 | contains
1541
1542		function FF_UsesDirectionalAtoms() result(doesit)
1543		logical :: doesit
1544	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1360	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1361	<	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1544	>	doesit = FF_uses_DirectionalAtoms
1545		end function FF_UsesDirectionalAtoms
1546
1547		function FF_RequiresPrepairCalc() result(doesit)
1548		logical :: doesit
1549	<	doesit = FF_uses_EAM
1549	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1550	>	.or. FF_uses_MEAM
1551		end function FF_RequiresPrepairCalc
1552
1369	–	function FF_RequiresPostpairCalc() result(doesit)
1370	–	logical :: doesit
1371	–	doesit = FF_uses_RF
1372	–	end function FF_RequiresPostpairCalc
1373	–
1553		#ifdef PROFILE
1554		function getforcetime() result(totalforcetime)
1555		real(kind=dp) :: totalforcetime

Diff Legend

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