ViewVC Help

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

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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2267 by tim, Fri Jul 29 17:34:06 2005 UTC vs.
Revision 2917 by chrisfen, Mon Jul 3 13:18:43 2006 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.25 2005-07-29 17:34:06 tim Exp $, $Date: 2005-07-29 17:34:06 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
48	>	!! @version $Id: doForces.F90,v 1.84 2006-07-03 13:18:43 chrisfen Exp $, $Date: 2006-07-03 13:18:43 $, $Name: not supported by cvs2svn $, $Revision: 1.84 $
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 72 | Line 72 | module doForces
72		PRIVATE
73
74		#define __FORTRAN90
75	<	#include "UseTheForce/fSwitchingFunction.h"
75	>	#include "UseTheForce/fCutoffPolicy.h"
76		#include "UseTheForce/DarkSide/fInteractionMap.h"
77	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78
79		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
80		INTEGER, PARAMETER:: PAIR_LOOP    = 2
81
81	–	logical, save :: haveRlist = .false.
82		logical, save :: haveNeighborList = .false.
83		logical, save :: haveSIMvariables = .false.
84		logical, save :: haveSaneForceField = .false.
85	<	logical, save :: haveInteractionMap = .false.
85	>	logical, save :: haveInteractionHash = .false.
86	>	logical, save :: haveGtypeCutoffMap = .false.
87	>	logical, save :: haveDefaultCutoffs = .false.
88	>	logical, save :: haveSkinThickness = .false.
89	>	logical, save :: haveElectrostaticSummationMethod = .false.
90	>	logical, save :: haveCutoffPolicy = .false.
91	>	logical, save :: VisitCutoffsAfterComputing = .false.
92	>	logical, save :: do_box_dipole = .false.
93
94		logical, save :: FF_uses_DirectionalAtoms
88	–	logical, save :: FF_uses_LennardJones
89	–	logical, save :: FF_uses_Electrostatics
90	–	logical, save :: FF_uses_Charges
95		logical, save :: FF_uses_Dipoles
92	–	logical, save :: FF_uses_Quadrupoles
93	–	logical, save :: FF_uses_Sticky
94	–	logical, save :: FF_uses_StickyPower
96		logical, save :: FF_uses_GayBerne
97		logical, save :: FF_uses_EAM
98	<	logical, save :: FF_uses_Shapes
99	<	logical, save :: FF_uses_FLARB
100	<	logical, save :: FF_uses_RF
98	>	logical, save :: FF_uses_SC
99	>	logical, save :: FF_uses_MEAM
100	>
101
102		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
103		logical, save :: SIM_uses_EAM
104	<	logical, save :: SIM_uses_Shapes
105	<	logical, save :: SIM_uses_FLARB
113	<	logical, save :: SIM_uses_RF
104	>	logical, save :: SIM_uses_SC
105	>	logical, save :: SIM_uses_MEAM
106		logical, save :: SIM_requires_postpair_calc
107		logical, save :: SIM_requires_prepair_calc
108		logical, save :: SIM_uses_PBC
117	–	logical, save :: SIM_uses_molecular_cutoffs
109
110	<	!!!GO AWAY---------
111	<	!!!!!real(kind=dp), save :: rlist, rlistsq
110	>	integer, save :: electrostaticSummationMethod
111	>	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
112
113	+	real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
114	+	real(kind=dp), save :: skinThickness
115	+	logical, save :: defaultDoShift
116	+
117		public :: init_FF
118	+	public :: setCutoffs
119	+	public :: cWasLame
120	+	public :: setElectrostaticMethod
121	+	public :: setBoxDipole
122	+	public :: getBoxDipole
123	+	public :: setCutoffPolicy
124	+	public :: setSkinThickness
125		public :: do_force_loop
124	–	!  public :: setRlistDF
125	–	!public :: addInteraction
126	–	!public :: setInteractionHash
127	–	!public :: getInteractionHash
128	–	public :: createInteractionMap
129	–	public :: createRcuts
126
127		#ifdef PROFILE
128		public :: getforcetime
#	Line 134 | Line 130 | module doForces
130		real :: forceTimeInitial, forceTimeFinal
131		integer :: nLoops
132		#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
133
134	<	type(Interaction), dimension(:,:),allocatable :: InteractionMap
135	<
134	>	!! Variables for cutoff mapping and interaction mapping
135	>	! Bit hash to determine pair-pair interactions.
136	>	integer, dimension(:,:), allocatable :: InteractionHash
137	>	real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
138	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
139	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
140
141	<
141	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
142	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
143	>
144	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
145	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
146	>	type ::gtypeCutoffs
147	>	real(kind=dp) :: rcut
148	>	real(kind=dp) :: rcutsq
149	>	real(kind=dp) :: rlistsq
150	>	end type gtypeCutoffs
151	>	type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
152	>
153	>	real(kind=dp), dimension(3) :: boxDipole
154	>
155		contains
156
157	<
151	<	subroutine createInteractionMap(status)
157	>	subroutine createInteractionHash()
158		integer :: nAtypes
153	–	integer, intent(out) :: status
159		integer :: i
160		integer :: j
161	<	integer :: ihash
157	<	real(kind=dp) :: myRcut
161	>	integer :: iHash
162		!! Test Types
163		logical :: i_is_LJ
164		logical :: i_is_Elect
#	Line 163 | Line 167 | contains
167		logical :: i_is_GB
168		logical :: i_is_EAM
169		logical :: i_is_Shape
170	+	logical :: i_is_SC
171	+	logical :: i_is_MEAM
172		logical :: j_is_LJ
173		logical :: j_is_Elect
174		logical :: j_is_Sticky
#	Line 170 | Line 176 | contains
176		logical :: j_is_GB
177		logical :: j_is_EAM
178		logical :: j_is_Shape
179	<
180	<	status = 0
181	<
179	>	logical :: j_is_SC
180	>	logical :: j_is_MEAM
181	>	real(kind=dp) :: myRcut
182	>
183		if (.not. associated(atypes)) then
184	<	call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
178	<	status = -1
184	>	call handleError("doForces", "atypes was not present before call of createInteractionHash!")
185		return
186		endif
187
188		nAtypes = getSize(atypes)
189
190		if (nAtypes == 0) then
191	<	status = -1
191	>	call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
192		return
193		end if
194
195	<	if (.not. allocated(InteractionMap)) then
196	<	allocate(InteractionMap(nAtypes,nAtypes))
195	>	if (.not. allocated(InteractionHash)) then
196	>	allocate(InteractionHash(nAtypes,nAtypes))
197	>	else
198	>	deallocate(InteractionHash)
199	>	allocate(InteractionHash(nAtypes,nAtypes))
200		endif
201	+
202	+	if (.not. allocated(atypeMaxCutoff)) then
203	+	allocate(atypeMaxCutoff(nAtypes))
204	+	else
205	+	deallocate(atypeMaxCutoff)
206	+	allocate(atypeMaxCutoff(nAtypes))
207	+	endif
208
209		do i = 1, nAtypes
210		call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
#	Line 198 | Line 214 | contains
214		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
215		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
216		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
217	+	call getElementProperty(atypes, i, "is_SC", i_is_SC)
218	+	call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
219
220		do j = i, nAtypes
221
#	Line 211 | Line 229 | contains
229		call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
230		call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
231		call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
232	+	call getElementProperty(atypes, j, "is_SC", j_is_SC)
233	+	call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
234
235		if (i_is_LJ .and. j_is_LJ) then
236		iHash = ior(iHash, LJ_PAIR)
#	Line 232 | Line 252 | contains
252		iHash = ior(iHash, EAM_PAIR)
253		endif
254
255	+	if (i_is_SC .and. j_is_SC) then
256	+	iHash = ior(iHash, SC_PAIR)
257	+	endif
258	+
259		if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
260		if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
261		if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
#	Line 241 | Line 265 | contains
265		if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
266
267
268	<	InteractionMap(i,j)%InteractionHash = iHash
269	<	InteractionMap(j,i)%InteractionHash = iHash
268	>	InteractionHash(i,j) = iHash
269	>	InteractionHash(j,i) = iHash
270
271		end do
272
273		end do
274
275	<	haveInteractionMap = .true.
276	<	end subroutine createInteractionMap
275	>	haveInteractionHash = .true.
276	>	end subroutine createInteractionHash
277
278	<	! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff.
255	<	subroutine createRcuts(defaultRList,stat)
256	<	real(kind=dp), intent(in), optional :: defaultRList
257	<	integer :: iMap
258	<	integer :: map_i,map_j
259	<	real(kind=dp) :: thisRCut = 0.0_dp
260	<	real(kind=dp) :: actualCutoff = 0.0_dp
261	<	integer, intent(out) :: stat
262	<	integer :: nAtypes
263	<	integer :: myStatus
278	>	subroutine createGtypeCutoffMap()
279
280	<	stat = 0
281	<	if (.not. haveInteractionMap) then
280	>	logical :: i_is_LJ
281	>	logical :: i_is_Elect
282	>	logical :: i_is_Sticky
283	>	logical :: i_is_StickyP
284	>	logical :: i_is_GB
285	>	logical :: i_is_EAM
286	>	logical :: i_is_Shape
287	>	logical :: i_is_SC
288	>	logical :: GtypeFound
289
290	<	call createInteractionMap(myStatus)
290	>	integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
291	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
292	>	integer :: nGroupsInRow
293	>	integer :: nGroupsInCol
294	>	integer :: nGroupTypesRow,nGroupTypesCol
295	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
296	>	real(kind=dp) :: biggestAtypeCutoff
297
298	<	if (myStatus .ne. 0) then
299	<	write(default_error, *) 'createInteractionMap failed in doForces!'
300	<	stat = -1
301	<	return
298	>	if (.not. haveInteractionHash) then
299	>	call createInteractionHash()
300	>	endif
301	>	#ifdef IS_MPI
302	>	nGroupsInRow = getNgroupsInRow(plan_group_row)
303	>	nGroupsInCol = getNgroupsInCol(plan_group_col)
304	>	#endif
305	>	nAtypes = getSize(atypes)
306	>	! Set all of the initial cutoffs to zero.
307	>	atypeMaxCutoff = 0.0_dp
308	>	do i = 1, nAtypes
309	>	if (SimHasAtype(i)) then
310	>	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
311	>	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
312	>	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
313	>	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
314	>	call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
315	>	call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
316	>	call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
317	>	call getElementProperty(atypes, i, "is_SC", i_is_SC)
318	>
319	>	if (haveDefaultCutoffs) then
320	>	atypeMaxCutoff(i) = defaultRcut
321	>	else
322	>	if (i_is_LJ) then
323	>	thisRcut = getSigma(i) * 2.5_dp
324	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325	>	endif
326	>	if (i_is_Elect) then
327	>	thisRcut = defaultRcut
328	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329	>	endif
330	>	if (i_is_Sticky) then
331	>	thisRcut = getStickyCut(i)
332	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333	>	endif
334	>	if (i_is_StickyP) then
335	>	thisRcut = getStickyPowerCut(i)
336	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337	>	endif
338	>	if (i_is_GB) then
339	>	thisRcut = getGayBerneCut(i)
340	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
341	>	endif
342	>	if (i_is_EAM) then
343	>	thisRcut = getEAMCut(i)
344	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
345	>	endif
346	>	if (i_is_Shape) then
347	>	thisRcut = getShapeCut(i)
348	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
349	>	endif
350	>	if (i_is_SC) then
351	>	thisRcut = getSCCut(i)
352	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
353	>	endif
354	>	endif
355	>
356	>	if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
357	>	biggestAtypeCutoff = atypeMaxCutoff(i)
358	>	endif
359	>
360		endif
361	+	enddo
362	+
363	+	istart = 1
364	+	jstart = 1
365	+	#ifdef IS_MPI
366	+	iend = nGroupsInRow
367	+	jend = nGroupsInCol
368	+	#else
369	+	iend = nGroups
370	+	jend = nGroups
371	+	#endif
372	+
373	+	!! allocate the groupToGtype and gtypeMaxCutoff here.
374	+	if(.not.allocated(groupToGtypeRow)) then
375	+	!  allocate(groupToGtype(iend))
376	+	allocate(groupToGtypeRow(iend))
377	+	else
378	+	deallocate(groupToGtypeRow)
379	+	allocate(groupToGtypeRow(iend))
380		endif
381	<
382	<
383	<	nAtypes = getSize(atypes)
384	<	!! If we pass a default rcut, set all atypes to that cutoff distance
385	<	if(present(defaultRList)) then
281	<	InteractionMap(:,:)%rList = defaultRList
282	<	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
283	<	haveRlist = .true.
284	<	return
381	>	if(.not.allocated(groupMaxCutoffRow)) then
382	>	allocate(groupMaxCutoffRow(iend))
383	>	else
384	>	deallocate(groupMaxCutoffRow)
385	>	allocate(groupMaxCutoffRow(iend))
386		end if
387
388	<	do map_i = 1,nAtypes
389	<	do map_j = map_i,nAtypes
390	<	iMap = InteractionMap(map_i, map_j)%InteractionHash
391	<
392	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
393	<	! thisRCut = getLJCutOff(map_i,map_j)
293	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
294	<	endif
295	<
296	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
297	<	! thisRCut = getElectrostaticCutOff(map_i,map_j)
298	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
299	<	endif
300	<
301	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
302	<	! thisRCut = getStickyCutOff(map_i,map_j)
303	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
304	<	endif
305	<
306	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
307	<	! thisRCut = getStickyPowerCutOff(map_i,map_j)
308	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
309	<	endif
310	<
311	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
312	<	! thisRCut = getGayberneCutOff(map_i,map_j)
313	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
314	<	endif
315	<
316	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
317	<	!              thisRCut = getGaybrneLJCutOff(map_i,map_j)
318	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
319	<	endif
320	<
321	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
322	<	!              thisRCut = getEAMCutOff(map_i,map_j)
323	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
324	<	endif
325	<
326	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
327	<	!              thisRCut = getShapeCutOff(map_i,map_j)
328	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
329	<	endif
330	<
331	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
332	<	!              thisRCut = getShapeLJCutOff(map_i,map_j)
333	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
334	<	endif
335	<	InteractionMap(map_i, map_j)%rList = actualCutoff
336	<	InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
337	<	end do
338	<	end do
339	<	haveRlist = .true.
340	<	end subroutine createRcuts
388	>	if(.not.allocated(gtypeMaxCutoffRow)) then
389	>	allocate(gtypeMaxCutoffRow(iend))
390	>	else
391	>	deallocate(gtypeMaxCutoffRow)
392	>	allocate(gtypeMaxCutoffRow(iend))
393	>	endif
394
395
396	<	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
397	<	!!$  subroutine setRlistDF( this_rlist )
398	<	!!$
399	<	!!$   real(kind=dp) :: this_rlist
400	<	!!$
401	<	!!$    rlist = this_rlist
402	<	!!$    rlistsq = rlist * rlist
403	<	!!$
351	<	!!$    haveRlist = .true.
352	<	!!$
353	<	!!$  end subroutine setRlistDF
396	>	#ifdef IS_MPI
397	>	! We only allocate new storage if we are in MPI because Ncol /= Nrow
398	>	if(.not.associated(groupToGtypeCol)) then
399	>	allocate(groupToGtypeCol(jend))
400	>	else
401	>	deallocate(groupToGtypeCol)
402	>	allocate(groupToGtypeCol(jend))
403	>	end if
404
405	<
406	<	subroutine setSimVariables()
407	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
408	<	SIM_uses_LennardJones = SimUsesLennardJones()
409	<	SIM_uses_Electrostatics = SimUsesElectrostatics()
410	<	SIM_uses_Charges = SimUsesCharges()
411	<	SIM_uses_Dipoles = SimUsesDipoles()
412	<	SIM_uses_Sticky = SimUsesSticky()
413	<	SIM_uses_StickyPower = SimUsesStickyPower()
414	<	SIM_uses_GayBerne = SimUsesGayBerne()
415	<	SIM_uses_EAM = SimUsesEAM()
416	<	SIM_uses_Shapes = SimUsesShapes()
367	<	SIM_uses_FLARB = SimUsesFLARB()
368	<	SIM_uses_RF = SimUsesRF()
369	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
370	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
371	<	SIM_uses_PBC = SimUsesPBC()
372	<
373	<	haveSIMvariables = .true.
405	>	if(.not.associated(groupMaxCutoffCol)) then
406	>	allocate(groupMaxCutoffCol(jend))
407	>	else
408	>	deallocate(groupMaxCutoffCol)
409	>	allocate(groupMaxCutoffCol(jend))
410	>	end if
411	>	if(.not.associated(gtypeMaxCutoffCol)) then
412	>	allocate(gtypeMaxCutoffCol(jend))
413	>	else
414	>	deallocate(gtypeMaxCutoffCol)
415	>	allocate(gtypeMaxCutoffCol(jend))
416	>	end if
417
418	<	return
419	<	end subroutine setSimVariables
418	>	groupMaxCutoffCol = 0.0_dp
419	>	gtypeMaxCutoffCol = 0.0_dp
420
421	<	subroutine doReadyCheck(error)
422	<	integer, intent(out) :: error
421	>	#endif
422	>	groupMaxCutoffRow = 0.0_dp
423	>	gtypeMaxCutoffRow = 0.0_dp
424
381	–	integer :: myStatus
425
426	<	error = 0
426	>	!! first we do a single loop over the cutoff groups to find the
427	>	!! largest cutoff for any atypes present in this group.  We also
428	>	!! create gtypes at this point.
429	>
430	>	tol = 1.0e-6_dp
431	>	nGroupTypesRow = 0
432	>	nGroupTypesCol = 0
433	>	do i = istart, iend
434	>	n_in_i = groupStartRow(i+1) - groupStartRow(i)
435	>	groupMaxCutoffRow(i) = 0.0_dp
436	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
437	>	atom1 = groupListRow(ia)
438	>	#ifdef IS_MPI
439	>	me_i = atid_row(atom1)
440	>	#else
441	>	me_i = atid(atom1)
442	>	#endif
443	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
444	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
445	>	endif
446	>	enddo
447	>	if (nGroupTypesRow.eq.0) then
448	>	nGroupTypesRow = nGroupTypesRow + 1
449	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
450	>	groupToGtypeRow(i) = nGroupTypesRow
451	>	else
452	>	GtypeFound = .false.
453	>	do g = 1, nGroupTypesRow
454	>	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
455	>	groupToGtypeRow(i) = g
456	>	GtypeFound = .true.
457	>	endif
458	>	enddo
459	>	if (.not.GtypeFound) then
460	>	nGroupTypesRow = nGroupTypesRow + 1
461	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
462	>	groupToGtypeRow(i) = nGroupTypesRow
463	>	endif
464	>	endif
465	>	enddo
466
467	<	if (.not. haveInteractionMap) then
467	>	#ifdef IS_MPI
468	>	do j = jstart, jend
469	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
470	>	groupMaxCutoffCol(j) = 0.0_dp
471	>	do ja = groupStartCol(j), groupStartCol(j+1)-1
472	>	atom1 = groupListCol(ja)
473
474	<	myStatus = 0
474	>	me_j = atid_col(atom1)
475
476	<	call createInteractionMap(myStatus)
476	>	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
477	>	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
478	>	endif
479	>	enddo
480
481	<	if (myStatus .ne. 0) then
482	<	write(default_error, *) 'createInteractionMap failed in doForces!'
483	<	error = -1
484	<	return
481	>	if (nGroupTypesCol.eq.0) then
482	>	nGroupTypesCol = nGroupTypesCol + 1
483	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
484	>	groupToGtypeCol(j) = nGroupTypesCol
485	>	else
486	>	GtypeFound = .false.
487	>	do g = 1, nGroupTypesCol
488	>	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
489	>	groupToGtypeCol(j) = g
490	>	GtypeFound = .true.
491	>	endif
492	>	enddo
493	>	if (.not.GtypeFound) then
494	>	nGroupTypesCol = nGroupTypesCol + 1
495	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
496	>	groupToGtypeCol(j) = nGroupTypesCol
497	>	endif
498		endif
499	+	enddo
500	+
501	+	#else
502	+	! Set pointers to information we just found
503	+	nGroupTypesCol = nGroupTypesRow
504	+	groupToGtypeCol => groupToGtypeRow
505	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
506	+	groupMaxCutoffCol => groupMaxCutoffRow
507	+	#endif
508	+
509	+	!! allocate the gtypeCutoffMap here.
510	+	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
511	+	!! then we do a double loop over all the group TYPES to find the cutoff
512	+	!! map between groups of two types
513	+	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
514	+
515	+	do i = 1, nGroupTypesRow
516	+	do j = 1, nGroupTypesCol
517	+
518	+	select case(cutoffPolicy)
519	+	case(TRADITIONAL_CUTOFF_POLICY)
520	+	thisRcut = tradRcut
521	+	case(MIX_CUTOFF_POLICY)
522	+	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
523	+	case(MAX_CUTOFF_POLICY)
524	+	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
525	+	case default
526	+	call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
527	+	return
528	+	end select
529	+	gtypeCutoffMap(i,j)%rcut = thisRcut
530	+
531	+	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
532	+
533	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
534	+
535	+	if (.not.haveSkinThickness) then
536	+	skinThickness = 1.0_dp
537	+	endif
538	+
539	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
540	+
541	+	! sanity check
542	+
543	+	if (haveDefaultCutoffs) then
544	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
545	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
546	+	endif
547	+	endif
548	+	enddo
549	+	enddo
550	+
551	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
552	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
553	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
554	+	#ifdef IS_MPI
555	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
556	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
557	+	#endif
558	+	groupMaxCutoffCol => null()
559	+	gtypeMaxCutoffCol => null()
560	+
561	+	haveGtypeCutoffMap = .true.
562	+	end subroutine createGtypeCutoffMap
563	+
564	+	subroutine setCutoffs(defRcut, defRsw)
565	+
566	+	real(kind=dp),intent(in) :: defRcut, defRsw
567	+	character(len = statusMsgSize) :: errMsg
568	+	integer :: localError
569	+
570	+	defaultRcut = defRcut
571	+	defaultRsw = defRsw
572	+
573	+	defaultDoShift = .false.
574	+	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
575	+
576	+	write(errMsg, *) &
577	+	'cutoffRadius and switchingRadius are set to the same', newline &
578	+	// tab, 'value.  OOPSE will use shifted ', newline &
579	+	// tab, 'potentials instead of switching functions.'
580	+
581	+	call handleInfo("setCutoffs", errMsg)
582	+
583	+	defaultDoShift = .true.
584	+
585	+	endif
586	+
587	+	localError = 0
588	+	call setLJDefaultCutoff( defaultRcut, defaultDoShift )
589	+	call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
590	+	call setCutoffEAM( defaultRcut )
591	+	call setCutoffSC( defaultRcut )
592	+	call set_switch(defaultRsw, defaultRcut)
593	+	call setHmatDangerousRcutValue(defaultRcut)
594	+
595	+	haveDefaultCutoffs = .true.
596	+	haveGtypeCutoffMap = .false.
597	+
598	+	end subroutine setCutoffs
599	+
600	+	subroutine cWasLame()
601	+
602	+	VisitCutoffsAfterComputing = .true.
603	+	return
604	+
605	+	end subroutine cWasLame
606	+
607	+	subroutine setCutoffPolicy(cutPolicy)
608	+
609	+	integer, intent(in) :: cutPolicy
610	+
611	+	cutoffPolicy = cutPolicy
612	+	haveCutoffPolicy = .true.
613	+	haveGtypeCutoffMap = .false.
614	+
615	+	end subroutine setCutoffPolicy
616	+
617	+	subroutine setBoxDipole()
618	+
619	+	do_box_dipole = .true.
620	+
621	+	end subroutine setBoxDipole
622	+
623	+	subroutine getBoxDipole( box_dipole )
624	+
625	+	real(kind=dp), intent(inout), dimension(3) :: box_dipole
626	+
627	+	box_dipole = boxDipole
628	+
629	+	end subroutine getBoxDipole
630	+
631	+	subroutine setElectrostaticMethod( thisESM )
632	+
633	+	integer, intent(in) :: thisESM
634	+
635	+	electrostaticSummationMethod = thisESM
636	+	haveElectrostaticSummationMethod = .true.
637	+
638	+	end subroutine setElectrostaticMethod
639	+
640	+	subroutine setSkinThickness( thisSkin )
641	+
642	+	real(kind=dp), intent(in) :: thisSkin
643	+
644	+	skinThickness = thisSkin
645	+	haveSkinThickness = .true.
646	+	haveGtypeCutoffMap = .false.
647	+
648	+	end subroutine setSkinThickness
649	+
650	+	subroutine setSimVariables()
651	+	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
652	+	SIM_uses_EAM = SimUsesEAM()
653	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
654	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
655	+	SIM_uses_PBC = SimUsesPBC()
656	+	SIM_uses_SC = SimUsesSC()
657	+
658	+	haveSIMvariables = .true.
659	+
660	+	return
661	+	end subroutine setSimVariables
662	+
663	+	subroutine doReadyCheck(error)
664	+	integer, intent(out) :: error
665	+	integer :: myStatus
666	+
667	+	error = 0
668	+
669	+	if (.not. haveInteractionHash) then
670	+	call createInteractionHash()
671		endif
672
673	<	if (.not. haveSIMvariables) then
674	<	call setSimVariables()
673	>	if (.not. haveGtypeCutoffMap) then
674	>	call createGtypeCutoffMap()
675		endif
676
677	<	if (.not. haveRlist) then
678	<	write(default_error, *) 'rList has not been set in doForces!'
679	<	error = -1
680	<	return
677	>	if (VisitCutoffsAfterComputing) then
678	>	call set_switch(largestRcut, largestRcut)
679	>	call setHmatDangerousRcutValue(largestRcut)
680	>	call setCutoffEAM(largestRcut)
681	>	call setCutoffSC(largestRcut)
682	>	VisitCutoffsAfterComputing = .false.
683		endif
684
685	+	if (.not. haveSIMvariables) then
686	+	call setSimVariables()
687	+	endif
688	+
689		if (.not. haveNeighborList) then
690		write(default_error, *) 'neighbor list has not been initialized in doForces!'
691		error = -1
692		return
693		end if
694	<
694	>
695		if (.not. haveSaneForceField) then
696		write(default_error, *) 'Force Field is not sane in doForces!'
697		error = -1
698		return
699		end if
700	<
700	>
701		#ifdef IS_MPI
702		if (.not. isMPISimSet()) then
703		write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
#	Line 428 | Line 709 | contains
709		end subroutine doReadyCheck
710
711
712	<	subroutine init_FF(use_RF_c, thisStat)
712	>	subroutine init_FF(thisStat)
713
433	–	logical, intent(in) :: use_RF_c
434	–
714		integer, intent(out) :: thisStat
715		integer :: my_status, nMatches
716		integer, pointer :: MatchList(:) => null()
438	–	real(kind=dp) :: rcut, rrf, rt, dielect
717
718		!! assume things are copacetic, unless they aren't
719		thisStat = 0
720
443	–	!! Fortran's version of a cast:
444	–	FF_uses_RF = use_RF_c
445	–
721		!! init_FF is called *after* all of the atom types have been
722		!! defined in atype_module using the new_atype subroutine.
723		!!
#	Line 450 | Line 725 | contains
725		!! interactions are used by the force field.
726
727		FF_uses_DirectionalAtoms = .false.
453	–	FF_uses_LennardJones = .false.
454	–	FF_uses_Electrostatics = .false.
455	–	FF_uses_Charges = .false.
728		FF_uses_Dipoles = .false.
457	–	FF_uses_Sticky = .false.
458	–	FF_uses_StickyPower = .false.
729		FF_uses_GayBerne = .false.
730		FF_uses_EAM = .false.
731	<	FF_uses_Shapes = .false.
462	<	FF_uses_FLARB = .false.
731	>	FF_uses_SC = .false.
732
733		call getMatchingElementList(atypes, "is_Directional", .true., &
734		nMatches, MatchList)
735		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
736
468	–	call getMatchingElementList(atypes, "is_LennardJones", .true., &
469	–	nMatches, MatchList)
470	–	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
471	–
472	–	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
473	–	nMatches, MatchList)
474	–	if (nMatches .gt. 0) then
475	–	FF_uses_Electrostatics = .true.
476	–	endif
477	–
478	–	call getMatchingElementList(atypes, "is_Charge", .true., &
479	–	nMatches, MatchList)
480	–	if (nMatches .gt. 0) then
481	–	FF_uses_Charges = .true.
482	–	FF_uses_Electrostatics = .true.
483	–	endif
484	–
737		call getMatchingElementList(atypes, "is_Dipole", .true., &
738		nMatches, MatchList)
739	<	if (nMatches .gt. 0) then
488	<	FF_uses_Dipoles = .true.
489	<	FF_uses_Electrostatics = .true.
490	<	FF_uses_DirectionalAtoms = .true.
491	<	endif
492	<
493	<	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
494	<	nMatches, MatchList)
495	<	if (nMatches .gt. 0) then
496	<	FF_uses_Quadrupoles = .true.
497	<	FF_uses_Electrostatics = .true.
498	<	FF_uses_DirectionalAtoms = .true.
499	<	endif
500	<
501	<	call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
502	<	MatchList)
503	<	if (nMatches .gt. 0) then
504	<	FF_uses_Sticky = .true.
505	<	FF_uses_DirectionalAtoms = .true.
506	<	endif
507	<
508	<	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
509	<	MatchList)
510	<	if (nMatches .gt. 0) then
511	<	FF_uses_StickyPower = .true.
512	<	FF_uses_DirectionalAtoms = .true.
513	<	endif
739	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
740
741		call getMatchingElementList(atypes, "is_GayBerne", .true., &
742		nMatches, MatchList)
743	<	if (nMatches .gt. 0) then
518	<	FF_uses_GayBerne = .true.
519	<	FF_uses_DirectionalAtoms = .true.
520	<	endif
743	>	if (nMatches .gt. 0) FF_uses_GayBerne = .true.
744
745		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
746		if (nMatches .gt. 0) FF_uses_EAM = .true.
747
748	<	call getMatchingElementList(atypes, "is_Shape", .true., &
749	<	nMatches, MatchList)
527	<	if (nMatches .gt. 0) then
528	<	FF_uses_Shapes = .true.
529	<	FF_uses_DirectionalAtoms = .true.
530	<	endif
748	>	call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
749	>	if (nMatches .gt. 0) FF_uses_SC = .true.
750
532	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
533	–	nMatches, MatchList)
534	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
751
536	–	!! Assume sanity (for the sake of argument)
752		haveSaneForceField = .true.
753
539	–	!! check to make sure the FF_uses_RF setting makes sense
540	–
541	–	if (FF_uses_dipoles) then
542	–	if (FF_uses_RF) then
543	–	dielect = getDielect()
544	–	call initialize_rf(dielect)
545	–	endif
546	–	else
547	–	if (FF_uses_RF) then
548	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
549	–	thisStat = -1
550	–	haveSaneForceField = .false.
551	–	return
552	–	endif
553	–	endif
554	–
555	–	!sticky module does not contain check_sticky_FF anymore
556	–	!if (FF_uses_sticky) then
557	–	!   call check_sticky_FF(my_status)
558	–	!   if (my_status /= 0) then
559	–	!      thisStat = -1
560	–	!      haveSaneForceField = .false.
561	–	!      return
562	–	!   end if
563	–	!endif
564	–
754		if (FF_uses_EAM) then
755		call init_EAM_FF(my_status)
756		if (my_status /= 0) then
#	Line 572 | Line 761 | contains
761		end if
762		endif
763
575	–	if (FF_uses_GayBerne) then
576	–	call check_gb_pair_FF(my_status)
577	–	if (my_status .ne. 0) then
578	–	thisStat = -1
579	–	haveSaneForceField = .false.
580	–	return
581	–	endif
582	–	endif
583	–
584	–	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
585	–	endif
586	–
764		if (.not. haveNeighborList) then
765		!! Create neighbor lists
766		call expandNeighborList(nLocal, my_status)
#	Line 617 | Line 794 | contains
794
795		!! Stress Tensor
796		real( kind = dp), dimension(9) :: tau
797	<	real ( kind = dp ) :: pot
797	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
798		logical ( kind = 2) :: do_pot_c, do_stress_c
799		logical :: do_pot
800		logical :: do_stress
801		logical :: in_switching_region
802		#ifdef IS_MPI
803	<	real( kind = DP ) :: pot_local
803	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
804		integer :: nAtomsInRow
805		integer :: nAtomsInCol
806		integer :: nprocs
#	Line 638 | Line 815 | contains
815		integer :: nlist
816		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
817		real( kind = DP ) :: sw, dswdr, swderiv, mf
818	+	real( kind = DP ) :: rVal
819		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
820		real(kind=dp) :: rfpot, mu_i, virial
821	+	real(kind=dp):: rCut
822		integer :: me_i, me_j, n_in_i, n_in_j
823		logical :: is_dp_i
824		integer :: neighborListSize
#	Line 647 | Line 826 | contains
826		integer :: localError
827		integer :: propPack_i, propPack_j
828		integer :: loopStart, loopEnd, loop
829	<	integer :: iMap
830	<	real(kind=dp) :: listSkin = 1.0
829	>	integer :: iHash
830	>	integer :: i1
831
832	+	!! the variables for the box dipole moment
833	+	#ifdef IS_MPI
834	+	integer :: pChgCount_local
835	+	integer :: nChgCount_local
836	+	real(kind=dp) :: pChg_local
837	+	real(kind=dp) :: nChg_local
838	+	real(kind=dp), dimension(3) :: pChgPos_local
839	+	real(kind=dp), dimension(3) :: nChgPos_local
840	+	real(kind=dp), dimension(3) :: dipVec_local
841	+	#endif
842	+	integer :: pChgCount
843	+	integer :: nChgCount
844	+	real(kind=dp) :: pChg
845	+	real(kind=dp) :: nChg
846	+	real(kind=dp) :: chg_value
847	+	real(kind=dp), dimension(3) :: pChgPos
848	+	real(kind=dp), dimension(3) :: nChgPos
849	+	real(kind=dp), dimension(3) :: dipVec
850	+	real(kind=dp), dimension(3) :: chgVec
851	+
852	+	!! initialize box dipole variables
853	+	if (do_box_dipole) then
854	+	#ifdef IS_MPI
855	+	pChg_local = 0.0_dp
856	+	nChg_local = 0.0_dp
857	+	pChgCount_local = 0
858	+	nChgCount_local = 0
859	+	do i=1, 3
860	+	pChgPos_local = 0.0_dp
861	+	nChgPos_local = 0.0_dp
862	+	dipVec_local = 0.0_dp
863	+	enddo
864	+	#endif
865	+	pChg = 0.0_dp
866	+	nChg = 0.0_dp
867	+	pChgCount = 0
868	+	nChgCount = 0
869	+	chg_value = 0.0_dp
870	+
871	+	do i=1, 3
872	+	pChgPos(i) = 0.0_dp
873	+	nChgPos(i) = 0.0_dp
874	+	dipVec(i) = 0.0_dp
875	+	chgVec(i) = 0.0_dp
876	+	boxDipole(i) = 0.0_dp
877	+	enddo
878	+	endif
879	+
880		!! initialize local variables
881
882		#ifdef IS_MPI
#	Line 712 | Line 939 | contains
939		! (but only on the first time through):
940		if (loop .eq. loopStart) then
941		#ifdef IS_MPI
942	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
942	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
943		update_nlist)
944		#else
945	<	call checkNeighborList(nGroups, q_group, listSkin, &
945	>	call checkNeighborList(nGroups, q_group, skinThickness, &
946		update_nlist)
947		#endif
948		endif
#	Line 739 | Line 966 | contains
966		#endif
967		outer: do i = istart, iend
968
742	–	#ifdef IS_MPI
743	–	me_i = atid_row(i)
744	–	#else
745	–	me_i = atid(i)
746	–	#endif
747	–
969		if (update_nlist) point(i) = nlist + 1
970
971		n_in_i = groupStartRow(i+1) - groupStartRow(i)
#	Line 779 | Line 1000 | contains
1000		me_j = atid(j)
1001		call get_interatomic_vector(q_group(:,i), &
1002		q_group(:,j), d_grp, rgrpsq)
1003	<	#endif
1003	>	#endif
1004
1005	<	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
1005	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
1006		if (update_nlist) then
1007		nlist = nlist + 1
1008
#	Line 801 | Line 1022 | contains
1022
1023		list(nlist) = j
1024		endif
1025	<
1026	<	if (loop .eq. PAIR_LOOP) then
806	<	vij = 0.0d0
807	<	fij(1:3) = 0.0d0
808	<	endif
1025	>
1026	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
1027
1028	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
1029	<	in_switching_region)
1030	<
1031	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
1032	<
1033	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
1034	<
1035	<	atom1 = groupListRow(ia)
1036	<
1037	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1038	<
1039	<	atom2 = groupListCol(jb)
1040	<
1041	<	if (skipThisPair(atom1, atom2)) cycle inner
1042	<
1043	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1044	<	d_atm(1:3) = d_grp(1:3)
1045	<	ratmsq = rgrpsq
1046	<	else
1028	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
1029	>	if (loop .eq. PAIR_LOOP) then
1030	>	vij = 0.0_dp
1031	>	fij(1) = 0.0_dp
1032	>	fij(2) = 0.0_dp
1033	>	fij(3) = 0.0_dp
1034	>	endif
1035	>
1036	>	call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
1037	>
1038	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
1039	>
1040	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
1041	>
1042	>	atom1 = groupListRow(ia)
1043	>
1044	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1045	>
1046	>	atom2 = groupListCol(jb)
1047	>
1048	>	if (skipThisPair(atom1, atom2))  cycle inner
1049	>
1050	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1051	>	d_atm(1) = d_grp(1)
1052	>	d_atm(2) = d_grp(2)
1053	>	d_atm(3) = d_grp(3)
1054	>	ratmsq = rgrpsq
1055	>	else
1056		#ifdef IS_MPI
1057	<	call get_interatomic_vector(q_Row(:,atom1), &
1058	<	q_Col(:,atom2), d_atm, ratmsq)
1057	>	call get_interatomic_vector(q_Row(:,atom1), &
1058	>	q_Col(:,atom2), d_atm, ratmsq)
1059		#else
1060	<	call get_interatomic_vector(q(:,atom1), &
1061	<	q(:,atom2), d_atm, ratmsq)
1060	>	call get_interatomic_vector(q(:,atom1), &
1061	>	q(:,atom2), d_atm, ratmsq)
1062		#endif
1063	<	endif
1064	<
1065	<	if (loop .eq. PREPAIR_LOOP) then
1063	>	endif
1064	>
1065	>	if (loop .eq. PREPAIR_LOOP) then
1066		#ifdef IS_MPI
1067	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1068	<	rgrpsq, d_grp, do_pot, do_stress, &
1069	<	eFrame, A, f, t, pot_local)
1067	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1068	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1069	>	eFrame, A, f, t, pot_local)
1070		#else
1071	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1072	<	rgrpsq, d_grp, do_pot, do_stress, &
1073	<	eFrame, A, f, t, pot)
1071	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1072	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1073	>	eFrame, A, f, t, pot)
1074		#endif
1075	<	else
1075	>	else
1076		#ifdef IS_MPI
1077	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1078	<	do_pot, &
1079	<	eFrame, A, f, t, pot_local, vpair, fpair)
1077	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1078	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1079	>	fpair, d_grp, rgrp, rCut)
1080		#else
1081	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1082	<	do_pot,  &
1083	<	eFrame, A, f, t, pot, vpair, fpair)
1081	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1082	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1083	>	d_grp, rgrp, rCut)
1084		#endif
1085	+	vij = vij + vpair
1086	+	fij(1) = fij(1) + fpair(1)
1087	+	fij(2) = fij(2) + fpair(2)
1088	+	fij(3) = fij(3) + fpair(3)
1089	+	endif
1090	+	enddo inner
1091	+	enddo
1092
1093	<	vij = vij + vpair
1094	<	fij(1:3) = fij(1:3) + fpair(1:3)
1095	<	endif
1096	<	enddo inner
1097	<	enddo
1098	<
1099	<	if (loop .eq. PAIR_LOOP) then
1100	<	if (in_switching_region) then
1101	<	swderiv = vij*dswdr/rgrp
1102	<	fij(1) = fij(1) + swderiv*d_grp(1)
869	<	fij(2) = fij(2) + swderiv*d_grp(2)
870	<	fij(3) = fij(3) + swderiv*d_grp(3)
871	<
872	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
873	<	atom1=groupListRow(ia)
874	<	mf = mfactRow(atom1)
1093	>	if (loop .eq. PAIR_LOOP) then
1094	>	if (in_switching_region) then
1095	>	swderiv = vij*dswdr/rgrp
1096	>	fij(1) = fij(1) + swderiv*d_grp(1)
1097	>	fij(2) = fij(2) + swderiv*d_grp(2)
1098	>	fij(3) = fij(3) + swderiv*d_grp(3)
1099	>
1100	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1101	>	atom1=groupListRow(ia)
1102	>	mf = mfactRow(atom1)
1103		#ifdef IS_MPI
1104	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1105	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1106	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1104	>	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1105	>	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1106	>	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1107		#else
1108	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1109	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1110	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1108	>	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1109	>	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1110	>	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1111		#endif
1112	<	enddo
1113	<
1114	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1115	<	atom2=groupListCol(jb)
1116	<	mf = mfactCol(atom2)
1112	>	enddo
1113	>
1114	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1115	>	atom2=groupListCol(jb)
1116	>	mf = mfactCol(atom2)
1117		#ifdef IS_MPI
1118	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1119	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1120	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1118	>	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1119	>	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1120	>	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1121		#else
1122	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1123	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1124	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1122	>	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1123	>	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1124	>	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1125		#endif
1126	<	enddo
1127	<	endif
1126	>	enddo
1127	>	endif
1128
1129	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1129	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1130	>	endif
1131		endif
1132	<	end if
1132	>	endif
1133		enddo
1134	+
1135		enddo outer
1136
1137		if (update_nlist) then
#	Line 961 | Line 1191 | contains
1191
1192		if (do_pot) then
1193		! scatter/gather pot_row into the members of my column
1194	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1195	<
1194	>	do i = 1,LR_POT_TYPES
1195	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1196	>	end do
1197		! scatter/gather pot_local into all other procs
1198		! add resultant to get total pot
1199		do i = 1, nlocal
1200	<	pot_local = pot_local + pot_Temp(i)
1200	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1201	>	+ pot_Temp(1:LR_POT_TYPES,i)
1202		enddo
1203
1204		pot_Temp = 0.0_DP
1205	<
1206	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1205	>	do i = 1,LR_POT_TYPES
1206	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1207	>	end do
1208		do i = 1, nlocal
1209	<	pot_local = pot_local + pot_Temp(i)
1209	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1210	>	+ pot_Temp(1:LR_POT_TYPES,i)
1211		enddo
1212
1213		endif
1214		#endif
1215
1216	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1216	>	if (SIM_requires_postpair_calc) then
1217	>	do i = 1, nlocal
1218	>
1219	>	! we loop only over the local atoms, so we don't need row and column
1220	>	! lookups for the types
1221	>
1222	>	me_i = atid(i)
1223	>
1224	>	! is the atom electrostatic?  See if it would have an
1225	>	! electrostatic interaction with itself
1226	>	iHash = InteractionHash(me_i,me_i)
1227
1228	<	if (FF_uses_RF .and. SIM_uses_RF) then
985	<
1228	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1229		#ifdef IS_MPI
1230	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1231	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
989	<	do i = 1,nlocal
990	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
991	<	end do
992	<	#endif
993	<
994	<	do i = 1, nLocal
995	<
996	<	rfpot = 0.0_DP
997	<	#ifdef IS_MPI
998	<	me_i = atid_row(i)
1230	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1231	>	t, do_pot)
1232		#else
1233	<	me_i = atid(i)
1233	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1234	>	t, do_pot)
1235		#endif
1236	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1236	>	endif
1237	>
1238	>
1239	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1240
1241	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1242	<
1243	<	mu_i = getDipoleMoment(me_i)
1244	<
1245	<	!! The reaction field needs to include a self contribution
1246	<	!! to the field:
1247	<	call accumulate_self_rf(i, mu_i, eFrame)
1248	<	!! Get the reaction field contribution to the
1249	<	!! potential and torques:
1250	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1241	>	! loop over the excludes to accumulate RF stuff we've
1242	>	! left out of the normal pair loop
1243	>
1244	>	do i1 = 1, nSkipsForAtom(i)
1245	>	j = skipsForAtom(i, i1)
1246	>
1247	>	! prevent overcounting of the skips
1248	>	if (i.lt.j) then
1249	>	call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1250	>	rVal = sqrt(ratmsq)
1251	>	call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1252		#ifdef IS_MPI
1253	<	pot_local = pot_local + rfpot
1253	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1254	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1255		#else
1256	<	pot = pot + rfpot
1256	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1257	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1258	>	#endif
1259	>	endif
1260	>	enddo
1261	>	endif
1262
1263	+	if (do_box_dipole) then
1264	+	#ifdef IS_MPI
1265	+	call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1266	+	nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1267	+	pChgCount_local, nChgCount_local)
1268	+	#else
1269	+	call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1270	+	pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1271		#endif
1272	<	endif
1273	<	enddo
1022	<	endif
1272	>	endif
1273	>	enddo
1274		endif
1275
1025	–
1276		#ifdef IS_MPI
1027	–
1277		if (do_pot) then
1278	<	pot = pot + pot_local
1279	<	!! we assume the c code will do the allreduce to get the total potential
1280	<	!! we could do it right here if we needed to...
1278	>	#ifdef SINGLE_PRECISION
1279	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1280	>	mpi_comm_world,mpi_err)
1281	>	#else
1282	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1283	>	mpi_sum, mpi_comm_world,mpi_err)
1284	>	#endif
1285		endif
1286	<
1286	>
1287		if (do_stress) then
1288	+	#ifdef SINGLE_PRECISION
1289	+	call mpi_allreduce(tau_Temp, tau, 9,mpi_real,mpi_sum, &
1290	+	mpi_comm_world,mpi_err)
1291	+	call mpi_allreduce(virial_Temp, virial,1,mpi_real,mpi_sum, &
1292	+	mpi_comm_world,mpi_err)
1293	+	#else
1294		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1295		mpi_comm_world,mpi_err)
1296		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1297		mpi_comm_world,mpi_err)
1298	+	#endif
1299		endif
1300	+
1301	+	if (do_box_dipole) then
1302
1303	+	#ifdef SINGLE_PRECISION
1304	+	call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1305	+	mpi_comm_world, mpi_err)
1306	+	call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1307	+	mpi_comm_world, mpi_err)
1308	+	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1309	+	mpi_comm_world, mpi_err)
1310	+	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1311	+	mpi_comm_world, mpi_err)
1312	+	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1313	+	mpi_comm_world, mpi_err)
1314	+	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1315	+	mpi_comm_world, mpi_err)
1316	+	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1317	+	mpi_comm_world, mpi_err)
1318		#else
1319	+	call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1320	+	mpi_comm_world, mpi_err)
1321	+	call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1322	+	mpi_comm_world, mpi_err)
1323	+	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1324	+	mpi_sum, mpi_comm_world, mpi_err)
1325	+	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1326	+	mpi_sum, mpi_comm_world, mpi_err)
1327	+	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1328	+	mpi_sum, mpi_comm_world, mpi_err)
1329	+	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1330	+	mpi_sum, mpi_comm_world, mpi_err)
1331	+	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1332	+	mpi_sum, mpi_comm_world, mpi_err)
1333	+	#endif
1334
1335	+	endif
1336	+
1337	+	#else
1338	+
1339		if (do_stress) then
1340		tau = tau_Temp
1341		virial = virial_Temp
1342		endif
1343	<
1343	>
1344		#endif
1345
1346	+	if (do_box_dipole) then
1347	+	! first load the accumulated dipole moment (if dipoles were present)
1348	+	boxDipole(1) = dipVec(1)
1349	+	boxDipole(2) = dipVec(2)
1350	+	boxDipole(3) = dipVec(3)
1351	+
1352	+	! now include the dipole moment due to charges
1353	+	! use the lesser of the positive and negative charge totals
1354	+	if (nChg .le. pChg) then
1355	+	chg_value = nChg
1356	+	else
1357	+	chg_value = pChg
1358	+	endif
1359	+
1360	+	! find the average positions
1361	+	if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1362	+	pChgPos = pChgPos / pChgCount
1363	+	nChgPos = nChgPos / nChgCount
1364	+	endif
1365	+
1366	+	! dipole is from the negative to the positive (physics notation)
1367	+	chgVec(1) = pChgPos(1) - nChgPos(1)
1368	+	chgVec(2) = pChgPos(2) - nChgPos(2)
1369	+	chgVec(3) = pChgPos(3) - nChgPos(3)
1370	+
1371	+	boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1372	+	boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1373	+	boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1374	+
1375	+	endif
1376	+
1377		end subroutine do_force_loop
1378
1379		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1380	<	eFrame, A, f, t, pot, vpair, fpair)
1380	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1381
1382	<	real( kind = dp ) :: pot, vpair, sw
1382	>	real( kind = dp ) :: vpair, sw
1383	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1384		real( kind = dp ), dimension(3) :: fpair
1385		real( kind = dp ), dimension(nLocal)   :: mfact
1386		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1063 | Line 1391 | contains
1391		logical, intent(inout) :: do_pot
1392		integer, intent(in) :: i, j
1393		real ( kind = dp ), intent(inout) :: rijsq
1394	<	real ( kind = dp )                :: r
1394	>	real ( kind = dp ), intent(inout) :: r_grp
1395		real ( kind = dp ), intent(inout) :: d(3)
1396	<	real ( kind = dp ) :: ebalance
1396	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1397	>	real ( kind = dp ), intent(inout) :: rCut
1398	>	real ( kind = dp ) :: r
1399	>	real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1400		integer :: me_i, me_j
1401	+	integer :: k
1402
1403	<	integer :: iMap
1403	>	integer :: iHash
1404
1405		r = sqrt(rijsq)
1406	<	vpair = 0.0d0
1407	<	fpair(1:3) = 0.0d0
1406	>
1407	>	vpair = 0.0_dp
1408	>	fpair(1:3) = 0.0_dp
1409
1410		#ifdef IS_MPI
1411		me_i = atid_row(i)
#	Line 1082 | Line 1415 | contains
1415		me_j = atid(j)
1416		#endif
1417
1418	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1419	<
1420	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1421	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1418	>	iHash = InteractionHash(me_i, me_j)
1419	>
1420	>	if ( iand(iHash, LJ_PAIR).ne.0 ) then
1421	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1422	>	pot(VDW_POT), f, do_pot)
1423		endif
1424	<
1425	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1426	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1427	<	pot, eFrame, f, t, do_pot)
1094	<
1095	<	if (FF_uses_RF .and. SIM_uses_RF) then
1096	<
1097	<	! CHECK ME (RF needs to know about all electrostatic types)
1098	<	call accumulate_rf(i, j, r, eFrame, sw)
1099	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1100	<	endif
1101	<
1424	>
1425	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1426	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1427	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1428		endif
1429	<
1430	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1429	>
1430	>	if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1431		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1432	<	pot, A, f, t, do_pot)
1432	>	pot(HB_POT), A, f, t, do_pot)
1433		endif
1434	<
1435	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1434	>
1435	>	if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1436		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1437	<	pot, A, f, t, do_pot)
1437	>	pot(HB_POT), A, f, t, do_pot)
1438		endif
1439	<
1440	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1439	>
1440	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1441		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1442	<	pot, A, f, t, do_pot)
1442	>	pot(VDW_POT), A, f, t, do_pot)
1443		endif
1444
1445	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1446	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1447	<	!           pot, A, f, t, do_pot)
1445	>	if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1446	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1447	>	pot(VDW_POT), A, f, t, do_pot)
1448		endif
1449	<
1450	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1451	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1452	<	do_pot)
1449	>
1450	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1451	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1452	>	pot(METALLIC_POT), f, do_pot)
1453		endif
1454	<
1455	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1454	>
1455	>	if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1456		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1457	<	pot, A, f, t, do_pot)
1457	>	pot(VDW_POT), A, f, t, do_pot)
1458		endif
1459	<
1460	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1459	>
1460	>	if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1461		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1462	<	pot, A, f, t, do_pot)
1462	>	pot(VDW_POT), A, f, t, do_pot)
1463		endif
1464	<
1464	>
1465	>	if ( iand(iHash, SC_PAIR).ne.0 ) then
1466	>	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1467	>	pot(METALLIC_POT), f, do_pot)
1468	>	endif
1469	>
1470		end subroutine do_pair
1471
1472	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1472	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1473		do_pot, do_stress, eFrame, A, f, t, pot)
1474
1475	<	real( kind = dp ) :: pot, sw
1475	>	real( kind = dp ) :: sw
1476	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1477		real( kind = dp ), dimension(9,nLocal) :: eFrame
1478		real (kind=dp), dimension(9,nLocal) :: A
1479		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1149 | Line 1481 | contains
1481
1482		logical, intent(inout) :: do_pot, do_stress
1483		integer, intent(in) :: i, j
1484	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1484	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1485		real ( kind = dp )                :: r, rc
1486		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1487
1488	<	integer :: me_i, me_j, iMap
1488	>	integer :: me_i, me_j, iHash
1489
1490	+	r = sqrt(rijsq)
1491	+
1492		#ifdef IS_MPI
1493		me_i = atid_row(i)
1494		me_j = atid_col(j)
#	Line 1163 | Line 1497 | contains
1497		me_j = atid(j)
1498		#endif
1499
1500	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1500	>	iHash = InteractionHash(me_i, me_j)
1501
1502	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1503	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1502	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1503	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1504		endif
1505	+
1506	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1507	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1508	+	endif
1509
1510		end subroutine do_prepair
1511
1512
1513		subroutine do_preforce(nlocal,pot)
1514		integer :: nlocal
1515	<	real( kind = dp ) :: pot
1515	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1516
1517		if (FF_uses_EAM .and. SIM_uses_EAM) then
1518	<	call calc_EAM_preforce_Frho(nlocal,pot)
1518	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1519		endif
1520	<
1521	<
1520	>	if (FF_uses_SC .and. SIM_uses_SC) then
1521	>	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1522	>	endif
1523		end subroutine do_preforce
1524
1525
#	Line 1192 | Line 1531 | contains
1531		real( kind = dp ) :: d(3), scaled(3)
1532		integer i
1533
1534	<	d(1:3) = q_j(1:3) - q_i(1:3)
1534	>	d(1) = q_j(1) - q_i(1)
1535	>	d(2) = q_j(2) - q_i(2)
1536	>	d(3) = q_j(3) - q_i(3)
1537
1538		! Wrap back into periodic box if necessary
1539		if ( SIM_uses_PBC ) then
1540
1541		if( .not.boxIsOrthorhombic ) then
1542		! calc the scaled coordinates.
1543	+	! scaled = matmul(HmatInv, d)
1544
1545	<	scaled = matmul(HmatInv, d)
1546	<
1545	>	scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1546	>	scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1547	>	scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1548	>
1549		! wrap the scaled coordinates
1550
1551	<	scaled = scaled  - anint(scaled)
1551	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1552	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1553	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1554
1209	–
1555		! calc the wrapped real coordinates from the wrapped scaled
1556		! coordinates
1557	+	! d = matmul(Hmat,scaled)
1558	+	d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1559	+	d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1560	+	d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1561
1213	–	d = matmul(Hmat,scaled)
1214	–
1562		else
1563		! calc the scaled coordinates.
1564
1565	<	do i = 1, 3
1566	<	scaled(i) = d(i) * HmatInv(i,i)
1565	>	scaled(1) = d(1) * HmatInv(1,1)
1566	>	scaled(2) = d(2) * HmatInv(2,2)
1567	>	scaled(3) = d(3) * HmatInv(3,3)
1568	>
1569	>	! wrap the scaled coordinates
1570	>
1571	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1572	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1573	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1574
1575	<	! wrap the scaled coordinates
1575	>	! calc the wrapped real coordinates from the wrapped scaled
1576	>	! coordinates
1577
1578	<	scaled(i) = scaled(i) - anint(scaled(i))
1578	>	d(1) = scaled(1)*Hmat(1,1)
1579	>	d(2) = scaled(2)*Hmat(2,2)
1580	>	d(3) = scaled(3)*Hmat(3,3)
1581
1225	–	! calc the wrapped real coordinates from the wrapped scaled
1226	–	! coordinates
1227	–
1228	–	d(i) = scaled(i)*Hmat(i,i)
1229	–	enddo
1582		endif
1583
1584		endif
1585
1586	<	r_sq = dot_product(d,d)
1586	>	r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1587
1588		end subroutine get_interatomic_vector
1589
#	Line 1263 | Line 1615 | contains
1615		pot_Col = 0.0_dp
1616		pot_Temp = 0.0_dp
1617
1266	–	rf_Row = 0.0_dp
1267	–	rf_Col = 0.0_dp
1268	–	rf_Temp = 0.0_dp
1269	–
1618		#endif
1619
1620		if (FF_uses_EAM .and. SIM_uses_EAM) then
1621		call clean_EAM()
1622		endif
1623
1276	–	rf = 0.0_dp
1624		tau_Temp = 0.0_dp
1625		virial_Temp = 0.0_dp
1626		end subroutine zero_work_arrays
#	Line 1362 | Line 1709 | contains
1709
1710		function FF_UsesDirectionalAtoms() result(doesit)
1711		logical :: doesit
1712	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1366	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1367	<	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1712	>	doesit = FF_uses_DirectionalAtoms
1713		end function FF_UsesDirectionalAtoms
1714
1715		function FF_RequiresPrepairCalc() result(doesit)
1716		logical :: doesit
1717	<	doesit = FF_uses_EAM
1717	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1718	>	.or. FF_uses_MEAM
1719		end function FF_RequiresPrepairCalc
1720
1375	–	function FF_RequiresPostpairCalc() result(doesit)
1376	–	logical :: doesit
1377	–	doesit = FF_uses_RF
1378	–	end function FF_RequiresPostpairCalc
1379	–
1721		#ifdef PROFILE
1722		function getforcetime() result(totalforcetime)
1723		real(kind=dp) :: totalforcetime

Diff Legend

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