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 2129 by chrisfen, Mon Mar 21 20:51:10 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.12 2005-03-21 20:51:06 chrisfen Exp $, $Date: 2005-03-21 20:51:06 $, $Name: not supported by cvs2svn $, $Revision: 1.12 $
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
80	–	logical, save :: haveRlist = .false.
84		logical, save :: haveNeighborList = .false.
85		logical, save :: haveSIMvariables = .false.
83	–	logical, save :: havePropertyMap = .false.
86		logical, save :: haveSaneForceField = .false.
87	<
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
87	–	logical, save :: FF_uses_LennardJones
88	–	logical, save :: FF_uses_Electrostatics
89	–	logical, save :: FF_uses_Charges
96		logical, save :: FF_uses_Dipoles
91	–	logical, save :: FF_uses_Quadrupoles
92	–	logical, save :: FF_uses_sticky
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
100	–	logical, save :: SIM_uses_LennardJones
101	–	logical, save :: SIM_uses_Electrostatics
102	–	logical, save :: SIM_uses_Charges
103	–	logical, save :: SIM_uses_Dipoles
104	–	logical, save :: SIM_uses_Quadrupoles
105	–	logical, save :: SIM_uses_Sticky
106	–	logical, save :: SIM_uses_GayBerne
104		logical, save :: SIM_uses_EAM
105	<	logical, save :: SIM_uses_Shapes
106	<	logical, save :: SIM_uses_FLARB
110	<	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
114	–	logical, save :: SIM_uses_molecular_cutoffs
110
111	<	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
120	–	public :: setRlistDF
125
126		#ifdef PROFILE
127		public :: getforcetime
#	Line 125 | Line 129 | module doForces
129		real :: forceTimeInitial, forceTimeFinal
130		integer :: nLoops
131		#endif
132	+
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	<	type :: Properties
141	<	logical :: is_Directional   = .false.
131	<	logical :: is_LennardJones  = .false.
132	<	logical :: is_Electrostatic = .false.
133	<	logical :: is_Charge        = .false.
134	<	logical :: is_Dipole        = .false.
135	<	logical :: is_Quadrupole    = .false.
136	<	logical :: is_Sticky        = .false.
137	<	logical :: is_GayBerne      = .false.
138	<	logical :: is_EAM           = .false.
139	<	logical :: is_Shape         = .false.
140	<	logical :: is_FLARB         = .false.
141	<	end type Properties
140	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
141	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
142
143	<	type(Properties), dimension(:),allocatable :: PropertyMap
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	<	subroutine setRlistDF( this_rlist )
148	<
149	<	real(kind=dp) :: this_rlist
150	<
151	<	rlist = this_rlist
152	<	rlistsq = rlist * rlist
153	<
154	<	haveRlist = .true.
155	<
156	<	end subroutine setRlistDF
157	<
158	<	subroutine createPropertyMap(status)
154	>	subroutine createInteractionHash()
155		integer :: nAtypes
160	–	integer :: status
156		integer :: i
157	<	logical :: thisProperty
158	<	real (kind=DP) :: thisDPproperty
157	>	integer :: j
158	>	integer :: iHash
159	>	!! Test Types
160	>	logical :: i_is_LJ
161	>	logical :: i_is_Elect
162	>	logical :: i_is_Sticky
163	>	logical :: i_is_StickyP
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
172	>	logical :: j_is_StickyP
173	>	logical :: j_is_GB
174	>	logical :: j_is_EAM
175	>	logical :: j_is_Shape
176	>	logical :: j_is_SC
177	>	logical :: j_is_MEAM
178	>	real(kind=dp) :: myRcut
179
180	<	status = 0
181	<
180	>	if (.not. associated(atypes)) then
181	>	call handleError("doForces", "atypes was not present before call of createInteractionHash!")
182	>	return
183	>	endif
184	>
185		nAtypes = getSize(atypes)
186	<
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(PropertyMap)) then
193	<	allocate(PropertyMap(nAtypes))
191	>
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_Directional", thisProperty)
208	<	PropertyMap(i)%is_Directional = thisProperty
207	>	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
208	>	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
209	>	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
210	>	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
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	<	call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
183	<	PropertyMap(i)%is_LennardJones = thisProperty
184	<
185	<	call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
186	<	PropertyMap(i)%is_Electrostatic = thisProperty
217	>	do j = i, nAtypes
218
219	<	call getElementProperty(atypes, i, "is_Charge", thisProperty)
220	<	PropertyMap(i)%is_Charge = thisProperty
190	<
191	<	call getElementProperty(atypes, i, "is_Dipole", thisProperty)
192	<	PropertyMap(i)%is_Dipole = thisProperty
219	>	iHash = 0
220	>	myRcut = 0.0_dp
221
222	<	call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
223	<	PropertyMap(i)%is_Quadrupole = thisProperty
222	>	call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
223	>	call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
224	>	call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
225	>	call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
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	<	call getElementProperty(atypes, i, "is_Sticky", thisProperty)
233	<	PropertyMap(i)%is_Sticky = thisProperty
232	>	if (i_is_LJ .and. j_is_LJ) then
233	>	iHash = ior(iHash, LJ_PAIR)
234	>	endif
235	>
236	>	if (i_is_Elect .and. j_is_Elect) then
237	>	iHash = ior(iHash, ELECTROSTATIC_PAIR)
238	>	endif
239	>
240	>	if (i_is_Sticky .and. j_is_Sticky) then
241	>	iHash = ior(iHash, STICKY_PAIR)
242	>	endif
243
244	<	call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
245	<	PropertyMap(i)%is_GayBerne = thisProperty
244	>	if (i_is_StickyP .and. j_is_StickyP) then
245	>	iHash = ior(iHash, STICKYPOWER_PAIR)
246	>	endif
247
248	<	call getElementProperty(atypes, i, "is_EAM", thisProperty)
249	<	PropertyMap(i)%is_EAM = thisProperty
248	>	if (i_is_EAM .and. j_is_EAM) then
249	>	iHash = ior(iHash, EAM_PAIR)
250	>	endif
251
252	<	call getElementProperty(atypes, i, "is_Shape", thisProperty)
253	<	PropertyMap(i)%is_Shape = thisProperty
252	>	if (i_is_SC .and. j_is_SC) then
253	>	iHash = ior(iHash, SC_PAIR)
254	>	endif
255
256	<	call getElementProperty(atypes, i, "is_FLARB", thisProperty)
257	<	PropertyMap(i)%is_FLARB = thisProperty
258	<	end do
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)
259
260	<	havePropertyMap = .true.
260	>	if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
261	>	if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
262	>	if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
263
215	–	end subroutine createPropertyMap
264
265	<	subroutine setSimVariables()
266	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
219	<	SIM_uses_LennardJones = SimUsesLennardJones()
220	<	SIM_uses_Electrostatics = SimUsesElectrostatics()
221	<	SIM_uses_Charges = SimUsesCharges()
222	<	SIM_uses_Dipoles = SimUsesDipoles()
223	<	SIM_uses_Sticky = SimUsesSticky()
224	<	SIM_uses_GayBerne = SimUsesGayBerne()
225	<	SIM_uses_EAM = SimUsesEAM()
226	<	SIM_uses_Shapes = SimUsesShapes()
227	<	SIM_uses_FLARB = SimUsesFLARB()
228	<	SIM_uses_RF = SimUsesRF()
229	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
230	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
231	<	SIM_uses_PBC = SimUsesPBC()
265	>	InteractionHash(i,j) = iHash
266	>	InteractionHash(j,i) = iHash
267
268	<	haveSIMvariables = .true.
268	>	end do
269
270	<	return
236	<	end subroutine setSimVariables
270	>	end do
271
272	<	subroutine doReadyCheck(error)
273	<	integer, intent(out) :: error
272	>	haveInteractionHash = .true.
273	>	end subroutine createInteractionHash
274
275	<	integer :: myStatus
275	>	subroutine createGtypeCutoffMap()
276
277	<	error = 0
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	>	! 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	>
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 (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
349	>	biggestAtypeCutoff = atypeMaxCutoff(i)
350	>	endif
351	>
352	>	endif
353	>	enddo
354
355	<	if (.not. havePropertyMap) then
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	<	myStatus = 0
380	>	if(.not.allocated(gtypeMaxCutoffRow)) then
381	>	allocate(gtypeMaxCutoffRow(iend))
382	>	else
383	>	deallocate(gtypeMaxCutoffRow)
384	>	allocate(gtypeMaxCutoffRow(iend))
385	>	endif
386
249	–	call createPropertyMap(myStatus)
387
388	<	if (myStatus .ne. 0) then
389	<	write(default_error, *) 'createPropertyMap failed in doForces!'
390	<	error = -1
391	<	return
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	>	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	+
646	+	integer :: myStatus
647	+
648	+	error = 0
649	+
650	+	if (.not. haveInteractionHash) then
651	+	call createInteractionHash()
652		endif
653
654	+	if (.not. haveGtypeCutoffMap) then
655	+	call createGtypeCutoffMap()
656	+	endif
657	+
658	+
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 286 | Line 692 | contains
692		#endif
693		return
694		end subroutine doReadyCheck
289	–
695
291	–	subroutine init_FF(use_RF_c, thisStat)
696
697	<	logical, intent(in) :: use_RF_c
697	>	subroutine init_FF(thisStat)
698
699		integer, intent(out) :: thisStat
700		integer :: my_status, nMatches
701		integer, pointer :: MatchList(:) => null()
298	–	real(kind=dp) :: rcut, rrf, rt, dielect
702
703		!! assume things are copacetic, unless they aren't
704		thisStat = 0
705
303	–	!! Fortran's version of a cast:
304	–	FF_uses_RF = use_RF_c
305	–
706		!! init_FF is called *after* all of the atom types have been
707		!! defined in atype_module using the new_atype subroutine.
708		!!
709		!! this will scan through the known atypes and figure out what
710		!! interactions are used by the force field.
711	<
711	>
712		FF_uses_DirectionalAtoms = .false.
313	–	FF_uses_LennardJones = .false.
314	–	FF_uses_Electrostatics = .false.
315	–	FF_uses_Charges = .false.
713		FF_uses_Dipoles = .false.
317	–	FF_uses_Sticky = .false.
714		FF_uses_GayBerne = .false.
715		FF_uses_EAM = .false.
716	<	FF_uses_Shapes = .false.
321	<	FF_uses_FLARB = .false.
322	<
716	>
717		call getMatchingElementList(atypes, "is_Directional", .true., &
718		nMatches, MatchList)
719		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
720
327	–	call getMatchingElementList(atypes, "is_LennardJones", .true., &
328	–	nMatches, MatchList)
329	–	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
330	–
331	–	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
332	–	nMatches, MatchList)
333	–	if (nMatches .gt. 0) then
334	–	FF_uses_Electrostatics = .true.
335	–	endif
336	–
337	–	call getMatchingElementList(atypes, "is_Charge", .true., &
338	–	nMatches, MatchList)
339	–	if (nMatches .gt. 0) then
340	–	FF_uses_Charges = .true.
341	–	FF_uses_Electrostatics = .true.
342	–	endif
343	–
721		call getMatchingElementList(atypes, "is_Dipole", .true., &
722		nMatches, MatchList)
723	<	if (nMatches .gt. 0) then
347	<	FF_uses_Dipoles = .true.
348	<	FF_uses_Electrostatics = .true.
349	<	FF_uses_DirectionalAtoms = .true.
350	<	endif
351	<
352	<	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
353	<	nMatches, MatchList)
354	<	if (nMatches .gt. 0) then
355	<	FF_uses_Quadrupoles = .true.
356	<	FF_uses_Electrostatics = .true.
357	<	FF_uses_DirectionalAtoms = .true.
358	<	endif
723	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
724
360	–	call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
361	–	MatchList)
362	–	if (nMatches .gt. 0) then
363	–	FF_uses_Sticky = .true.
364	–	FF_uses_DirectionalAtoms = .true.
365	–	endif
366	–
725		call getMatchingElementList(atypes, "is_GayBerne", .true., &
726		nMatches, MatchList)
727	<	if (nMatches .gt. 0) then
728	<	FF_uses_GayBerne = .true.
371	<	FF_uses_DirectionalAtoms = .true.
372	<	endif
373	<
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.
376	–
377	–	call getMatchingElementList(atypes, "is_Shape", .true., &
378	–	nMatches, MatchList)
379	–	if (nMatches .gt. 0) then
380	–	FF_uses_Shapes = .true.
381	–	FF_uses_DirectionalAtoms = .true.
382	–	endif
731
384	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
385	–	nMatches, MatchList)
386	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
732
388	–	!! Assume sanity (for the sake of argument)
733		haveSaneForceField = .true.
390	–
391	–	!! check to make sure the FF_uses_RF setting makes sense
392	–
393	–	if (FF_uses_dipoles) then
394	–	if (FF_uses_RF) then
395	–	dielect = getDielect()
396	–	call initialize_rf(dielect)
397	–	endif
398	–	else
399	–	if (FF_uses_RF) then
400	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
401	–	thisStat = -1
402	–	haveSaneForceField = .false.
403	–	return
404	–	endif
405	–	endif
734
407	–	!sticky module does not contain check_sticky_FF anymore
408	–	!if (FF_uses_sticky) then
409	–	!   call check_sticky_FF(my_status)
410	–	!   if (my_status /= 0) then
411	–	!      thisStat = -1
412	–	!      haveSaneForceField = .false.
413	–	!      return
414	–	!   end if
415	–	!endif
416	–
735		if (FF_uses_EAM) then
736	<	call init_EAM_FF(my_status)
736	>	call init_EAM_FF(my_status)
737		if (my_status /= 0) then
738		write(default_error, *) "init_EAM_FF returned a bad status"
739		thisStat = -1
#	Line 424 | Line 742 | contains
742		end if
743		endif
744
427	–	if (FF_uses_GayBerne) then
428	–	call check_gb_pair_FF(my_status)
429	–	if (my_status .ne. 0) then
430	–	thisStat = -1
431	–	haveSaneForceField = .false.
432	–	return
433	–	endif
434	–	endif
435	–
436	–	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
437	–	endif
438	–
745		if (.not. haveNeighborList) then
746		!! Create neighbor lists
747		call expandNeighborList(nLocal, my_status)
#	Line 445 | Line 751 | contains
751		return
752		endif
753		haveNeighborList = .true.
754	<	endif
755	<
754	>	endif
755	>
756		end subroutine init_FF
451	–
757
758	+
759		!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
760		!------------------------------------------------------------->
761		subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, &
#	Line 469 | 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 490 | 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 499 | Line 807 | contains
807		integer :: localError
808		integer :: propPack_i, propPack_j
809		integer :: loopStart, loopEnd, loop
810	+	integer :: iHash
811	+	integer :: i1
812	+
813
503	–	real(kind=dp) :: listSkin = 1.0
504	–
814		!! initialize local variables
815	<
815	>
816		#ifdef IS_MPI
817		pot_local = 0.0_dp
818		nAtomsInRow   = getNatomsInRow(plan_atom_row)
#	Line 513 | Line 822 | contains
822		#else
823		natoms = nlocal
824		#endif
825	<
825	>
826		call doReadyCheck(localError)
827		if ( localError .ne. 0 ) then
828		call handleError("do_force_loop", "Not Initialized")
#	Line 521 | Line 830 | contains
830		return
831		end if
832		call zero_work_arrays()
833	<
833	>
834		do_pot = do_pot_c
835		do_stress = do_stress_c
836	<
836	>
837		! Gather all information needed by all force loops:
838	<
838	>
839		#ifdef IS_MPI
840	<
840	>
841		call gather(q, q_Row, plan_atom_row_3d)
842		call gather(q, q_Col, plan_atom_col_3d)
843
844		call gather(q_group, q_group_Row, plan_group_row_3d)
845		call gather(q_group, q_group_Col, plan_group_col_3d)
846	<
846	>
847		if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
848		call gather(eFrame, eFrame_Row, plan_atom_row_rotation)
849		call gather(eFrame, eFrame_Col, plan_atom_col_rotation)
850	<
850	>
851		call gather(A, A_Row, plan_atom_row_rotation)
852		call gather(A, A_Col, plan_atom_col_rotation)
853		endif
854	<
854	>
855		#endif
856	<
856	>
857		!! Begin force loop timing:
858		#ifdef PROFILE
859		call cpu_time(forceTimeInitial)
860		nloops = nloops + 1
861		#endif
862	<
862	>
863		loopEnd = PAIR_LOOP
864		if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
865		loopStart = PREPAIR_LOOP
#	Line 564 | 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, &
877	<	update_nlist)
876	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
877	>	update_nlist)
878		#else
879	<	call checkNeighborList(nGroups, q_group, listSkin, &
880	<	update_nlist)
879	>	call checkNeighborList(nGroups, q_group, skinThickness, &
880	>	update_nlist)
881		#endif
882		endif
883	<
883	>
884		if (update_nlist) then
885		!! save current configuration and construct neighbor list
886		#ifdef IS_MPI
#	Line 582 | Line 891 | contains
891		neighborListSize = size(list)
892		nlist = 0
893		endif
894	<
894	>
895		istart = 1
896		#ifdef IS_MPI
897		iend = nGroupsInRow
#	Line 592 | Line 901 | contains
901		outer: do i = istart, iend
902
903		if (update_nlist) point(i) = nlist + 1
904	<
904	>
905		n_in_i = groupStartRow(i+1) - groupStartRow(i)
906	<
906	>
907		if (update_nlist) then
908		#ifdef IS_MPI
909		jstart = 1
#	Line 609 | Line 918 | contains
918		! make sure group i has neighbors
919		if (jstart .gt. jend) cycle outer
920		endif
921	<
921	>
922		do jnab = jstart, jend
923		if (update_nlist) then
924		j = jnab
#	Line 618 | Line 927 | contains
927		endif
928
929		#ifdef IS_MPI
930	+	me_j = atid_col(j)
931		call get_interatomic_vector(q_group_Row(:,i), &
932		q_group_Col(:,j), d_grp, rgrpsq)
933		#else
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 < rlistsq) then
939	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
940		if (update_nlist) then
941		nlist = nlist + 1
942	<
942	>
943		if (nlist > neighborListSize) then
944		#ifdef IS_MPI
945		call expandNeighborList(nGroupsInRow, listerror)
#	Line 642 | Line 953 | contains
953		end if
954		neighborListSize = size(list)
955		endif
956	<
956	>
957		list(nlist) = j
647	–	endif
648	–
649	–	if (loop .eq. PAIR_LOOP) then
650	–	vij = 0.0d0
651	–	fij(1:3) = 0.0d0
958		endif
959	+
960	+
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
962	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
963	>
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	<	atom1 = groupListRow(ia)
970	>	call get_switch(rgrpsq, sw, dswdr, rgrp, &
971	>	group_switch, in_switching_region)
972
973	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
973	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
974	>
975	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
976
977	<	atom2 = groupListCol(jb)
977	>	atom1 = groupListRow(ia)
978
979	<	if (skipThisPair(atom1, atom2)) cycle inner
980	<
981	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
982	<	d_atm(1:3) = d_grp(1:3)
983	<	ratmsq = rgrpsq
984	<	else
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)
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)
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
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	<
1019	<	vij = vij + vpair
1020	<	fij(1:3) = fij(1:3) + fpair(1:3)
1021	<	endif
1022	<	enddo inner
1023	<	enddo
1024	<
1025	<	if (loop .eq. PAIR_LOOP) then
1026	<	if (in_switching_region) then
1027	<	swderiv = vij*dswdr/rgrp
1028	<	fij(1) = fij(1) + swderiv*d_grp(1)
1029	<	fij(2) = fij(2) + swderiv*d_grp(2)
1030	<	fij(3) = fij(3) + swderiv*d_grp(3)
1031	<
1032	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
1033	<	atom1=groupListRow(ia)
718	<	mf = mfactRow(atom1)
1018	>	vij = vij + vpair
1019	>	fij(1:3) = fij(1:3) + fpair(1:3)
1020	>	endif
1021	>	enddo inner
1022	>	enddo
1023	>
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
1057	>	enddo
1058	>	endif
1059	>
1060	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1061		endif
744	–
745	–	if (do_stress) call add_stress_tensor(d_grp, fij)
1062		endif
1063	<	end if
1063	>	endif
1064		enddo
1065	+
1066		enddo outer
1067	<
1067	>
1068		if (update_nlist) then
1069		#ifdef IS_MPI
1070		point(nGroupsInRow + 1) = nlist + 1
#	Line 761 | Line 1078 | contains
1078		update_nlist = .false.
1079		endif
1080		endif
1081	<
1081	>
1082		if (loop .eq. PREPAIR_LOOP) then
1083		call do_preforce(nlocal, pot)
1084		endif
1085	<
1085	>
1086		enddo
1087	<
1087	>
1088		!! Do timing
1089		#ifdef PROFILE
1090		call cpu_time(forceTimeFinal)
1091		forceTime = forceTime + forceTimeFinal - forceTimeInitial
1092		#endif
1093	<
1093	>
1094		#ifdef IS_MPI
1095		!!distribute forces
1096	<
1096	>
1097		f_temp = 0.0_dp
1098		call scatter(f_Row,f_temp,plan_atom_row_3d)
1099		do i = 1,nlocal
1100		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
1101		end do
1102	<
1102	>
1103		f_temp = 0.0_dp
1104		call scatter(f_Col,f_temp,plan_atom_col_3d)
1105		do i = 1,nlocal
1106		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
1107		end do
1108	<
1108	>
1109		if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
1110		t_temp = 0.0_dp
1111		call scatter(t_Row,t_temp,plan_atom_row_3d)
#	Line 797 | Line 1114 | contains
1114		end do
1115		t_temp = 0.0_dp
1116		call scatter(t_Col,t_temp,plan_atom_col_3d)
1117	<
1117	>
1118		do i = 1,nlocal
1119		t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
1120		end do
1121		endif
1122	<
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	<
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	<
1143	>
1144		endif
1145		#endif
1146	<
1147	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1148	<
828	<	if (FF_uses_RF .and. SIM_uses_RF) then
1146	>
1147	>	if (SIM_requires_postpair_calc) then
1148	>	do i = 1, nlocal
1149
1150	<	#ifdef IS_MPI
1151	<	call scatter(rf_Row,rf,plan_atom_row_3d)
832	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
833	<	do i = 1,nlocal
834	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
835	<	end do
836	<	#endif
1150	>	! we loop only over the local atoms, so we don't need row and column
1151	>	! lookups for the types
1152
1153	<	do i = 1, nLocal
1154	<
1155	<	rfpot = 0.0_DP
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 ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1160		#ifdef IS_MPI
1161	<	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	+	endif
1168	+
1169	+
1170	+	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1171
1172	<	if (PropertyMap(me_i)%is_Dipole) then
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	<	mu_i = getDipoleMoment(me_i)
1179	<
1180	<	!! The reaction field needs to include a self contribution
1181	<	!! to the field:
1182	<	call accumulate_self_rf(i, mu_i, eFrame)
1183	<	!! Get the reaction field contribution to the
1184	<	!! potential and torques:
856	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
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
868	–
1198		#ifdef IS_MPI
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
874	<	!! 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
1205		if (do_stress) then
#	Line 889 | Line 1217 | contains
1217		endif
1218
1219		#endif
1220	<
1220	>
1221		end subroutine do_force_loop
1222	<
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 906 | 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 ), 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 :: iHash
1246	+
1247		r = sqrt(rijsq)
1248		vpair = 0.0d0
1249		fpair(1:3) = 0.0d0
#	Line 922 | Line 1256 | contains
1256		me_j = atid(j)
1257		#endif
1258
1259	<	!    write(*,*) i, j, me_i, me_j
1259	>	iHash = InteractionHash(me_i, me_j)
1260
1261	<	if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
1262	<
1263	<	if ( PropertyMap(me_i)%is_LennardJones .and. &
930	<	PropertyMap(me_j)%is_LennardJones ) then
931	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
932	<	endif
933	<
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 (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
1267	<
1268	<	if (PropertyMap(me_i)%is_Electrostatic .and. &
939	<	PropertyMap(me_j)%is_Electrostatic) then
940	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
941	<	pot, eFrame, f, t, do_pot)
942	<	endif
943	<
944	<	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
945	<	if ( PropertyMap(me_i)%is_Dipole .and. &
946	<	PropertyMap(me_j)%is_Dipole) then
947	<	if (FF_uses_RF .and. SIM_uses_RF) then
948	<	call accumulate_rf(i, j, r, eFrame, sw)
949	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
950	<	endif
951	<	endif
952	<	endif
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	<
1272	<	if (FF_uses_Sticky .and. SIM_uses_sticky) then
1273	<
958	<	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
959	<	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
960	<	pot, A, f, t, do_pot)
961	<	endif
962	<
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(HB_POT), A, f, t, do_pot)
1274		endif
1275	<
1276	<
1277	<	if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1278	<
968	<	if ( PropertyMap(me_i)%is_GayBerne .and. &
969	<	PropertyMap(me_j)%is_GayBerne) then
970	<	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
971	<	pot, A, f, t, do_pot)
972	<	endif
973	<
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(HB_POT), A, f, t, do_pot)
1279		endif
1280
1281	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1282	<
1283	<	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
979	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
980	<	do_pot)
981	<	endif
982	<
1281	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1282	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1283	>	pot(VDW_POT), A, f, t, do_pot)
1284		endif
1285	+
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(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(iHash, SHAPE_PAIR).ne.0 ) then
1297	+	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1298	+	pot(VDW_POT), A, f, t, do_pot)
1299	+	endif
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(VDW_POT), A, f, t, do_pot)
1304	+	endif
1305
1306	<
1307	<	!    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
1308	<
988	<	if (FF_uses_Shapes .and. SIM_uses_Shapes) then
989	<	if ( PropertyMap(me_i)%is_Shape .and. &
990	<	PropertyMap(me_j)%is_Shape ) then
991	<	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
992	<	pot, A, f, t, do_pot)
993	<	endif
994	<
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
1319	<	real( kind = dp ), dimension(9,nLocal) :: eFrame
1320	<	real (kind=dp), dimension(9,nLocal) :: A
1321	<	real (kind=dp), dimension(3,nLocal) :: f
1322	<	real (kind=dp), dimension(3,nLocal) :: t
1323	<
1008	<	logical, intent(inout) :: do_pot, do_stress
1009	<	integer, intent(in) :: i, j
1010	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1011	<	real ( kind = dp )                :: r, rc
1012	<	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1013	<
1014	<	logical :: is_EAM_i, is_EAM_j
1015	<
1016	<	integer :: me_i, me_j
1017	<
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
1323	>	real (kind=dp), dimension(3,nLocal) :: t
1324
1325	+	logical, intent(inout) :: do_pot, do_stress
1326	+	integer, intent(in) :: i, j
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, iHash
1332	+
1333		r = sqrt(rijsq)
1020	–	if (SIM_uses_molecular_cutoffs) then
1021	–	rc = sqrt(rcijsq)
1022	–	else
1023	–	rc = r
1024	–	endif
1025	–
1334
1335		#ifdef IS_MPI
1336	<	me_i = atid_row(i)
1337	<	me_j = atid_col(j)
1336	>	me_i = atid_row(i)
1337	>	me_j = atid_col(j)
1338		#else
1339	<	me_i = atid(i)
1340	<	me_j = atid(j)
1339	>	me_i = atid(i)
1340	>	me_j = atid(j)
1341		#endif
1342	<
1343	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1344	<
1345	<	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1346	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1347	<
1348	<	endif
1349	<
1350	<	end subroutine do_prepair
1351	<
1352	<
1353	<	subroutine do_preforce(nlocal,pot)
1354	<	integer :: nlocal
1355	<	real( kind = dp ) :: pot
1356	<
1357	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1358	<	call calc_EAM_preforce_Frho(nlocal,pot)
1359	<	endif
1360	<
1361	<
1362	<	end subroutine do_preforce
1363	<
1364	<
1365	<	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1366	<
1367	<	real (kind = dp), dimension(3) :: q_i
1368	<	real (kind = dp), dimension(3) :: q_j
1369	<	real ( kind = dp ), intent(out) :: r_sq
1370	<	real( kind = dp ) :: d(3), scaled(3)
1371	<	integer i
1372	<
1373	<	d(1:3) = q_j(1:3) - q_i(1:3)
1374	<
1375	<	! Wrap back into periodic box if necessary
1376	<	if ( SIM_uses_PBC ) then
1377	<
1378	<	if( .not.boxIsOrthorhombic ) then
1379	<	! calc the scaled coordinates.
1380	<
1381	<	scaled = matmul(HmatInv, d)
1382	<
1383	<	! wrap the scaled coordinates
1384	<
1385	<	scaled = scaled  - anint(scaled)
1386	<
1387	<
1388	<	! calc the wrapped real coordinates from the wrapped scaled
1389	<	! coordinates
1390	<
1391	<	d = matmul(Hmat,scaled)
1392	<
1393	<	else
1394	<	! calc the scaled coordinates.
1395	<
1396	<	do i = 1, 3
1397	<	scaled(i) = d(i) * HmatInv(i,i)
1398	<
1399	<	! wrap the scaled coordinates
1400	<
1401	<	scaled(i) = scaled(i) - anint(scaled(i))
1402	<
1403	<	! calc the wrapped real coordinates from the wrapped scaled
1404	<	! coordinates
1405	<
1406	<	d(i) = scaled(i)*Hmat(i,i)
1407	<	enddo
1408	<	endif
1409	<
1410	<	endif
1411	<
1412	<	r_sq = dot_product(d,d)
1413	<
1414	<	end subroutine get_interatomic_vector
1415	<
1416	<	subroutine zero_work_arrays()
1417	<
1342	>
1343	>	iHash = InteractionHash(me_i, me_j)
1344	>
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 ),dimension(LR_POT_TYPES) :: pot
1359	>
1360	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
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
1369	>
1370	>
1371	>	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1372	>
1373	>	real (kind = dp), dimension(3) :: q_i
1374	>	real (kind = dp), dimension(3) :: q_j
1375	>	real ( kind = dp ), intent(out) :: r_sq
1376	>	real( kind = dp ) :: d(3), scaled(3)
1377	>	integer i
1378	>
1379	>	d(1:3) = q_j(1:3) - q_i(1:3)
1380	>
1381	>	! Wrap back into periodic box if necessary
1382	>	if ( SIM_uses_PBC ) then
1383	>
1384	>	if( .not.boxIsOrthorhombic ) then
1385	>	! calc the scaled coordinates.
1386	>
1387	>	scaled = matmul(HmatInv, d)
1388	>
1389	>	! wrap the scaled coordinates
1390	>
1391	>	scaled = scaled  - anint(scaled)
1392	>
1393	>
1394	>	! calc the wrapped real coordinates from the wrapped scaled
1395	>	! coordinates
1396	>
1397	>	d = matmul(Hmat,scaled)
1398	>
1399	>	else
1400	>	! calc the scaled coordinates.
1401	>
1402	>	do i = 1, 3
1403	>	scaled(i) = d(i) * HmatInv(i,i)
1404	>
1405	>	! wrap the scaled coordinates
1406	>
1407	>	scaled(i) = scaled(i) - anint(scaled(i))
1408	>
1409	>	! calc the wrapped real coordinates from the wrapped scaled
1410	>	! coordinates
1411	>
1412	>	d(i) = scaled(i)*Hmat(i,i)
1413	>	enddo
1414	>	endif
1415	>
1416	>	endif
1417	>
1418	>	r_sq = dot_product(d,d)
1419	>
1420	>	end subroutine get_interatomic_vector
1421	>
1422	>	subroutine zero_work_arrays()
1423	>
1424		#ifdef IS_MPI
1111	–
1112	–	q_Row = 0.0_dp
1113	–	q_Col = 0.0_dp
1425
1426	<	q_group_Row = 0.0_dp
1427	<	q_group_Col = 0.0_dp
1428	<
1429	<	eFrame_Row = 0.0_dp
1430	<	eFrame_Col = 0.0_dp
1431	<
1432	<	A_Row = 0.0_dp
1433	<	A_Col = 0.0_dp
1434	<
1435	<	f_Row = 0.0_dp
1436	<	f_Col = 0.0_dp
1437	<	f_Temp = 0.0_dp
1438	<
1439	<	t_Row = 0.0_dp
1440	<	t_Col = 0.0_dp
1441	<	t_Temp = 0.0_dp
1442	<
1443	<	pot_Row = 0.0_dp
1444	<	pot_Col = 0.0_dp
1445	<	pot_Temp = 0.0_dp
1446	<
1447	<	rf_Row = 0.0_dp
1448	<	rf_Col = 0.0_dp
1449	<	rf_Temp = 0.0_dp
1139	<
1426	>	q_Row = 0.0_dp
1427	>	q_Col = 0.0_dp
1428	>
1429	>	q_group_Row = 0.0_dp
1430	>	q_group_Col = 0.0_dp
1431	>
1432	>	eFrame_Row = 0.0_dp
1433	>	eFrame_Col = 0.0_dp
1434	>
1435	>	A_Row = 0.0_dp
1436	>	A_Col = 0.0_dp
1437	>
1438	>	f_Row = 0.0_dp
1439	>	f_Col = 0.0_dp
1440	>	f_Temp = 0.0_dp
1441	>
1442	>	t_Row = 0.0_dp
1443	>	t_Col = 0.0_dp
1444	>	t_Temp = 0.0_dp
1445	>
1446	>	pot_Row = 0.0_dp
1447	>	pot_Col = 0.0_dp
1448	>	pot_Temp = 0.0_dp
1449	>
1450		#endif
1451	<
1452	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1453	<	call clean_EAM()
1454	<	endif
1455	<
1456	<	rf = 0.0_dp
1457	<	tau_Temp = 0.0_dp
1458	<	virial_Temp = 0.0_dp
1459	<	end subroutine zero_work_arrays
1460	<
1461	<	function skipThisPair(atom1, atom2) result(skip_it)
1462	<	integer, intent(in) :: atom1
1463	<	integer, intent(in), optional :: atom2
1464	<	logical :: skip_it
1465	<	integer :: unique_id_1, unique_id_2
1466	<	integer :: me_i,me_j
1467	<	integer :: i
1468	<
1469	<	skip_it = .false.
1470	<
1471	<	!! there are a number of reasons to skip a pair or a particle
1472	<	!! mostly we do this to exclude atoms who are involved in short
1473	<	!! range interactions (bonds, bends, torsions), but we also need
1474	<	!! to exclude some overcounted interactions that result from
1475	<	!! the parallel decomposition
1166	<
1451	>
1452	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1453	>	call clean_EAM()
1454	>	endif
1455	>
1456	>	tau_Temp = 0.0_dp
1457	>	virial_Temp = 0.0_dp
1458	>	end subroutine zero_work_arrays
1459	>
1460	>	function skipThisPair(atom1, atom2) result(skip_it)
1461	>	integer, intent(in) :: atom1
1462	>	integer, intent(in), optional :: atom2
1463	>	logical :: skip_it
1464	>	integer :: unique_id_1, unique_id_2
1465	>	integer :: me_i,me_j
1466	>	integer :: i
1467	>
1468	>	skip_it = .false.
1469	>
1470	>	!! there are a number of reasons to skip a pair or a particle
1471	>	!! mostly we do this to exclude atoms who are involved in short
1472	>	!! range interactions (bonds, bends, torsions), but we also need
1473	>	!! to exclude some overcounted interactions that result from
1474	>	!! the parallel decomposition
1475	>
1476		#ifdef IS_MPI
1477	<	!! in MPI, we have to look up the unique IDs for each atom
1478	<	unique_id_1 = AtomRowToGlobal(atom1)
1477	>	!! in MPI, we have to look up the unique IDs for each atom
1478	>	unique_id_1 = AtomRowToGlobal(atom1)
1479		#else
1480	<	!! in the normal loop, the atom numbers are unique
1481	<	unique_id_1 = atom1
1480	>	!! in the normal loop, the atom numbers are unique
1481	>	unique_id_1 = atom1
1482		#endif
1483	<
1484	<	!! We were called with only one atom, so just check the global exclude
1485	<	!! list for this atom
1486	<	if (.not. present(atom2)) then
1487	<	do i = 1, nExcludes_global
1488	<	if (excludesGlobal(i) == unique_id_1) then
1489	<	skip_it = .true.
1490	<	return
1491	<	end if
1492	<	end do
1493	<	return
1494	<	end if
1495	<
1483	>
1484	>	!! We were called with only one atom, so just check the global exclude
1485	>	!! list for this atom
1486	>	if (.not. present(atom2)) then
1487	>	do i = 1, nExcludes_global
1488	>	if (excludesGlobal(i) == unique_id_1) then
1489	>	skip_it = .true.
1490	>	return
1491	>	end if
1492	>	end do
1493	>	return
1494	>	end if
1495	>
1496		#ifdef IS_MPI
1497	<	unique_id_2 = AtomColToGlobal(atom2)
1497	>	unique_id_2 = AtomColToGlobal(atom2)
1498		#else
1499	<	unique_id_2 = atom2
1499	>	unique_id_2 = atom2
1500		#endif
1501	<
1501	>
1502		#ifdef IS_MPI
1503	<	!! this situation should only arise in MPI simulations
1504	<	if (unique_id_1 == unique_id_2) then
1505	<	skip_it = .true.
1506	<	return
1507	<	end if
1508	<
1509	<	!! this prevents us from doing the pair on multiple processors
1510	<	if (unique_id_1 < unique_id_2) then
1511	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1512	<	skip_it = .true.
1513	<	return
1514	<	endif
1515	<	else
1516	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1517	<	skip_it = .true.
1518	<	return
1519	<	endif
1520	<	endif
1503	>	!! this situation should only arise in MPI simulations
1504	>	if (unique_id_1 == unique_id_2) then
1505	>	skip_it = .true.
1506	>	return
1507	>	end if
1508	>
1509	>	!! this prevents us from doing the pair on multiple processors
1510	>	if (unique_id_1 < unique_id_2) then
1511	>	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1512	>	skip_it = .true.
1513	>	return
1514	>	endif
1515	>	else
1516	>	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1517	>	skip_it = .true.
1518	>	return
1519	>	endif
1520	>	endif
1521		#endif
1522	<
1523	<	!! the rest of these situations can happen in all simulations:
1524	<	do i = 1, nExcludes_global
1525	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1526	<	(excludesGlobal(i) == unique_id_2)) then
1527	<	skip_it = .true.
1528	<	return
1529	<	endif
1530	<	enddo
1531	<
1532	<	do i = 1, nSkipsForAtom(atom1)
1533	<	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1534	<	skip_it = .true.
1535	<	return
1536	<	endif
1537	<	end do
1538	<
1539	<	return
1540	<	end function skipThisPair
1541	<
1542	<	function FF_UsesDirectionalAtoms() result(doesit)
1543	<	logical :: doesit
1544	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1545	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1546	<	FF_uses_GayBerne .or. FF_uses_Shapes
1547	<	end function FF_UsesDirectionalAtoms
1548	<
1549	<	function FF_RequiresPrepairCalc() result(doesit)
1550	<	logical :: doesit
1551	<	doesit = FF_uses_EAM
1552	<	end function FF_RequiresPrepairCalc
1244	<
1245	<	function FF_RequiresPostpairCalc() result(doesit)
1246	<	logical :: doesit
1247	<	doesit = FF_uses_RF
1248	<	end function FF_RequiresPostpairCalc
1249	<
1522	>
1523	>	!! the rest of these situations can happen in all simulations:
1524	>	do i = 1, nExcludes_global
1525	>	if ((excludesGlobal(i) == unique_id_1) .or. &
1526	>	(excludesGlobal(i) == unique_id_2)) then
1527	>	skip_it = .true.
1528	>	return
1529	>	endif
1530	>	enddo
1531	>
1532	>	do i = 1, nSkipsForAtom(atom1)
1533	>	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1534	>	skip_it = .true.
1535	>	return
1536	>	endif
1537	>	end do
1538	>
1539	>	return
1540	>	end function skipThisPair
1541	>
1542	>	function FF_UsesDirectionalAtoms() result(doesit)
1543	>	logical :: doesit
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 .or. FF_uses_SC &
1550	>	.or. FF_uses_MEAM
1551	>	end function FF_RequiresPrepairCalc
1552	>
1553		#ifdef PROFILE
1554	<	function getforcetime() result(totalforcetime)
1555	<	real(kind=dp) :: totalforcetime
1556	<	totalforcetime = forcetime
1557	<	end function getforcetime
1554	>	function getforcetime() result(totalforcetime)
1555	>	real(kind=dp) :: totalforcetime
1556	>	totalforcetime = forcetime
1557	>	end function getforcetime
1558		#endif
1256	–
1257	–	!! This cleans componets of force arrays belonging only to fortran
1559
1560	<	subroutine add_stress_tensor(dpair, fpair)
1561	<
1562	<	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1563	<
1564	<	! because the d vector is the rj - ri vector, and
1565	<	! because fx, fy, fz are the force on atom i, we need a
1566	<	! negative sign here:
1567	<
1568	<	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1569	<	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1570	<	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1571	<	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1572	<	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1573	<	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1574	<	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1575	<	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1576	<	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1577	<
1578	<	virial_Temp = virial_Temp + &
1579	<	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1580	<
1581	<	end subroutine add_stress_tensor
1582	<
1560	>	!! This cleans componets of force arrays belonging only to fortran
1561	>
1562	>	subroutine add_stress_tensor(dpair, fpair)
1563	>
1564	>	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1565	>
1566	>	! because the d vector is the rj - ri vector, and
1567	>	! because fx, fy, fz are the force on atom i, we need a
1568	>	! negative sign here:
1569	>
1570	>	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1571	>	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1572	>	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1573	>	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1574	>	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1575	>	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1576	>	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1577	>	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1578	>	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1579	>
1580	>	virial_Temp = virial_Temp + &
1581	>	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1582	>
1583	>	end subroutine add_stress_tensor
1584	>
1585		end module doForces

Diff Legend

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