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 2722 by gezelter, Thu Apr 20 18:24:24 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.79 2006-04-20 18:24:24 gezelter Exp $, $Date: 2006-04-20 18:24:24 $, $Name: not supported by cvs2svn $, $Revision: 1.79 $
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	+	use interpolation
68		#ifdef IS_MPI
69		use mpiSimulation
70		#endif
#	Line 72 | Line 73 | module doForces
73		PRIVATE
74
75		#define __FORTRAN90
76	<	#include "UseTheForce/fSwitchingFunction.h"
76	>	#include "UseTheForce/fCutoffPolicy.h"
77		#include "UseTheForce/DarkSide/fInteractionMap.h"
78	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79
80		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81		INTEGER, PARAMETER:: PAIR_LOOP    = 2
82	+	INTEGER, PARAMETER:: np = 500
83
81	–	logical, save :: haveRlist = .false.
84		logical, save :: haveNeighborList = .false.
85		logical, save :: haveSIMvariables = .false.
86		logical, save :: haveSaneForceField = .false.
87	<	logical, save :: haveInteractionMap = .false.
87	>	logical, save :: haveInteractionHash = .false.
88	>	logical, save :: haveGtypeCutoffMap = .false.
89	>	logical, save :: haveDefaultCutoffs = .false.
90	>	logical, save :: haveSkinThickness = .false.
91	>	logical, save :: haveElectrostaticSummationMethod = .false.
92	>	logical, save :: haveCutoffPolicy = .false.
93	>	logical, save :: VisitCutoffsAfterComputing = .false.
94	>	logical, save :: haveSplineSqrt = .false.
95
96		logical, save :: FF_uses_DirectionalAtoms
88	–	logical, save :: FF_uses_LennardJones
89	–	logical, save :: FF_uses_Electrostatics
90	–	logical, save :: FF_uses_Charges
97		logical, save :: FF_uses_Dipoles
92	–	logical, save :: FF_uses_Quadrupoles
93	–	logical, save :: FF_uses_Sticky
94	–	logical, save :: FF_uses_StickyPower
98		logical, save :: FF_uses_GayBerne
99		logical, save :: FF_uses_EAM
100	<	logical, save :: FF_uses_Shapes
101	<	logical, save :: FF_uses_FLARB
102	<	logical, save :: FF_uses_RF
100	>	logical, save :: FF_uses_SC
101	>	logical, save :: FF_uses_MEAM
102	>
103
104		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
105		logical, save :: SIM_uses_EAM
106	<	logical, save :: SIM_uses_Shapes
107	<	logical, save :: SIM_uses_FLARB
113	<	logical, save :: SIM_uses_RF
106	>	logical, save :: SIM_uses_SC
107	>	logical, save :: SIM_uses_MEAM
108		logical, save :: SIM_requires_postpair_calc
109		logical, save :: SIM_requires_prepair_calc
110		logical, save :: SIM_uses_PBC
117	–	logical, save :: SIM_uses_molecular_cutoffs
111
112	<	!!!GO AWAY---------
113	<	!!!!!real(kind=dp), save :: rlist, rlistsq
112	>	integer, save :: electrostaticSummationMethod
113	>	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
114
115	+	real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
116	+	real(kind=dp), save :: skinThickness
117	+	logical, save :: defaultDoShift
118	+
119		public :: init_FF
120	+	public :: setCutoffs
121	+	public :: cWasLame
122	+	public :: setElectrostaticMethod
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	+	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	+	! variables for the spline of the sqrt
154	+	type(cubicSpline), save :: splineSqrt
155	+	logical, save :: useSpline = .true.
156
157	+
158		contains
159
160	<
151	<	subroutine createInteractionMap(status)
160	>	subroutine createInteractionHash()
161		integer :: nAtypes
153	–	integer, intent(out) :: status
162		integer :: i
163		integer :: j
164	<	integer :: ihash
157	<	real(kind=dp) :: myRcut
164	>	integer :: iHash
165		!! Test Types
166		logical :: i_is_LJ
167		logical :: i_is_Elect
#	Line 163 | Line 170 | contains
170		logical :: i_is_GB
171		logical :: i_is_EAM
172		logical :: i_is_Shape
173	+	logical :: i_is_SC
174	+	logical :: i_is_MEAM
175		logical :: j_is_LJ
176		logical :: j_is_Elect
177		logical :: j_is_Sticky
#	Line 170 | Line 179 | contains
179		logical :: j_is_GB
180		logical :: j_is_EAM
181		logical :: j_is_Shape
182	<
183	<	status = 0
184	<
182	>	logical :: j_is_SC
183	>	logical :: j_is_MEAM
184	>	real(kind=dp) :: myRcut
185	>
186		if (.not. associated(atypes)) then
187	<	call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
178	<	status = -1
187	>	call handleError("doForces", "atypes was not present before call of createInteractionHash!")
188		return
189		endif
190
191		nAtypes = getSize(atypes)
192
193		if (nAtypes == 0) then
194	<	status = -1
194	>	call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
195		return
196		end if
197
198	<	if (.not. allocated(InteractionMap)) then
199	<	allocate(InteractionMap(nAtypes,nAtypes))
198	>	if (.not. allocated(InteractionHash)) then
199	>	allocate(InteractionHash(nAtypes,nAtypes))
200	>	else
201	>	deallocate(InteractionHash)
202	>	allocate(InteractionHash(nAtypes,nAtypes))
203		endif
204	+
205	+	if (.not. allocated(atypeMaxCutoff)) then
206	+	allocate(atypeMaxCutoff(nAtypes))
207	+	else
208	+	deallocate(atypeMaxCutoff)
209	+	allocate(atypeMaxCutoff(nAtypes))
210	+	endif
211
212		do i = 1, nAtypes
213		call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
#	Line 198 | Line 217 | contains
217		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
218		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
219		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
220	+	call getElementProperty(atypes, i, "is_SC", i_is_SC)
221	+	call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
222
223		do j = i, nAtypes
224
#	Line 211 | Line 232 | contains
232		call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
233		call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
234		call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
235	+	call getElementProperty(atypes, j, "is_SC", j_is_SC)
236	+	call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
237
238		if (i_is_LJ .and. j_is_LJ) then
239		iHash = ior(iHash, LJ_PAIR)
#	Line 232 | Line 255 | contains
255		iHash = ior(iHash, EAM_PAIR)
256		endif
257
258	+	if (i_is_SC .and. j_is_SC) then
259	+	iHash = ior(iHash, SC_PAIR)
260	+	endif
261	+
262		if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
263		if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
264		if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
#	Line 241 | Line 268 | contains
268		if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
269
270
271	<	InteractionMap(i,j)%InteractionHash = iHash
272	<	InteractionMap(j,i)%InteractionHash = iHash
271	>	InteractionHash(i,j) = iHash
272	>	InteractionHash(j,i) = iHash
273
274		end do
275
276		end do
277
278	<	haveInteractionMap = .true.
279	<	end subroutine createInteractionMap
278	>	haveInteractionHash = .true.
279	>	end subroutine createInteractionHash
280
281	<	! 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
281	>	subroutine createGtypeCutoffMap()
282
283	<	stat = 0
284	<	if (.not. haveInteractionMap) then
283	>	logical :: i_is_LJ
284	>	logical :: i_is_Elect
285	>	logical :: i_is_Sticky
286	>	logical :: i_is_StickyP
287	>	logical :: i_is_GB
288	>	logical :: i_is_EAM
289	>	logical :: i_is_Shape
290	>	logical :: i_is_SC
291	>	logical :: GtypeFound
292
293	<	call createInteractionMap(myStatus)
293	>	integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
294	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
295	>	integer :: nGroupsInRow
296	>	integer :: nGroupsInCol
297	>	integer :: nGroupTypesRow,nGroupTypesCol
298	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
299	>	real(kind=dp) :: biggestAtypeCutoff
300
301	<	if (myStatus .ne. 0) then
302	<	write(default_error, *) 'createInteractionMap failed in doForces!'
272	<	stat = -1
273	<	return
274	<	endif
301	>	if (.not. haveInteractionHash) then
302	>	call createInteractionHash()
303		endif
304	<
305	<
304	>	#ifdef IS_MPI
305	>	nGroupsInRow = getNgroupsInRow(plan_group_row)
306	>	nGroupsInCol = getNgroupsInCol(plan_group_col)
307	>	#endif
308		nAtypes = getSize(atypes)
309	<	!! If we pass a default rcut, set all atypes to that cutoff distance
310	<	if(present(defaultRList)) then
311	<	InteractionMap(:,:)%rList = defaultRList
312	<	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
313	<	haveRlist = .true.
314	<	return
315	<	end if
316	<
317	<	do map_i = 1,nAtypes
318	<	do map_j = map_i,nAtypes
319	<	iMap = InteractionMap(map_i, map_j)%InteractionHash
320	<
321	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
322	<	! thisRCut = getLJCutOff(map_i,map_j)
323	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
309	>	! Set all of the initial cutoffs to zero.
310	>	atypeMaxCutoff = 0.0_dp
311	>	do i = 1, nAtypes
312	>	if (SimHasAtype(i)) then
313	>	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
314	>	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
315	>	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
316	>	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
317	>	call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
318	>	call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
319	>	call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
320	>	call getElementProperty(atypes, i, "is_SC", i_is_SC)
321	>
322	>	if (haveDefaultCutoffs) then
323	>	atypeMaxCutoff(i) = defaultRcut
324	>	else
325	>	if (i_is_LJ) then
326	>	thisRcut = getSigma(i) * 2.5_dp
327	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
328	>	endif
329	>	if (i_is_Elect) then
330	>	thisRcut = defaultRcut
331	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
332	>	endif
333	>	if (i_is_Sticky) then
334	>	thisRcut = getStickyCut(i)
335	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
336	>	endif
337	>	if (i_is_StickyP) then
338	>	thisRcut = getStickyPowerCut(i)
339	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
340	>	endif
341	>	if (i_is_GB) then
342	>	thisRcut = getGayBerneCut(i)
343	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
344	>	endif
345	>	if (i_is_EAM) then
346	>	thisRcut = getEAMCut(i)
347	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
348	>	endif
349	>	if (i_is_Shape) then
350	>	thisRcut = getShapeCut(i)
351	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
352	>	endif
353	>	if (i_is_SC) then
354	>	thisRcut = getSCCut(i)
355	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
356	>	endif
357		endif
358	<
359	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
360	<	! thisRCut = getElectrostaticCutOff(map_i,map_j)
298	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
358	>
359	>	if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
360	>	biggestAtypeCutoff = atypeMaxCutoff(i)
361		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
362
363	+	endif
364	+	enddo
365	+
366	+	istart = 1
367	+	jstart = 1
368	+	#ifdef IS_MPI
369	+	iend = nGroupsInRow
370	+	jend = nGroupsInCol
371	+	#else
372	+	iend = nGroups
373	+	jend = nGroups
374	+	#endif
375	+
376	+	!! allocate the groupToGtype and gtypeMaxCutoff here.
377	+	if(.not.allocated(groupToGtypeRow)) then
378	+	!  allocate(groupToGtype(iend))
379	+	allocate(groupToGtypeRow(iend))
380	+	else
381	+	deallocate(groupToGtypeRow)
382	+	allocate(groupToGtypeRow(iend))
383	+	endif
384	+	if(.not.allocated(groupMaxCutoffRow)) then
385	+	allocate(groupMaxCutoffRow(iend))
386	+	else
387	+	deallocate(groupMaxCutoffRow)
388	+	allocate(groupMaxCutoffRow(iend))
389	+	end if
390
391	<	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
392	<	!!$  subroutine setRlistDF( this_rlist )
393	<	!!$
394	<	!!$   real(kind=dp) :: this_rlist
395	<	!!$
396	<	!!$    rlist = this_rlist
349	<	!!$    rlistsq = rlist * rlist
350	<	!!$
351	<	!!$    haveRlist = .true.
352	<	!!$
353	<	!!$  end subroutine setRlistDF
391	>	if(.not.allocated(gtypeMaxCutoffRow)) then
392	>	allocate(gtypeMaxCutoffRow(iend))
393	>	else
394	>	deallocate(gtypeMaxCutoffRow)
395	>	allocate(gtypeMaxCutoffRow(iend))
396	>	endif
397
398
399	<	subroutine setSimVariables()
400	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
401	<	SIM_uses_LennardJones = SimUsesLennardJones()
402	<	SIM_uses_Electrostatics = SimUsesElectrostatics()
403	<	SIM_uses_Charges = SimUsesCharges()
404	<	SIM_uses_Dipoles = SimUsesDipoles()
405	<	SIM_uses_Sticky = SimUsesSticky()
406	<	SIM_uses_StickyPower = SimUsesStickyPower()
364	<	SIM_uses_GayBerne = SimUsesGayBerne()
365	<	SIM_uses_EAM = SimUsesEAM()
366	<	SIM_uses_Shapes = SimUsesShapes()
367	<	SIM_uses_FLARB = SimUsesFLARB()
368	<	SIM_uses_RF = SimUsesRF()
369	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
370	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
371	<	SIM_uses_PBC = SimUsesPBC()
399	>	#ifdef IS_MPI
400	>	! We only allocate new storage if we are in MPI because Ncol /= Nrow
401	>	if(.not.associated(groupToGtypeCol)) then
402	>	allocate(groupToGtypeCol(jend))
403	>	else
404	>	deallocate(groupToGtypeCol)
405	>	allocate(groupToGtypeCol(jend))
406	>	end if
407
408	<	haveSIMvariables = .true.
408	>	if(.not.associated(groupMaxCutoffCol)) then
409	>	allocate(groupMaxCutoffCol(jend))
410	>	else
411	>	deallocate(groupMaxCutoffCol)
412	>	allocate(groupMaxCutoffCol(jend))
413	>	end if
414	>	if(.not.associated(gtypeMaxCutoffCol)) then
415	>	allocate(gtypeMaxCutoffCol(jend))
416	>	else
417	>	deallocate(gtypeMaxCutoffCol)
418	>	allocate(gtypeMaxCutoffCol(jend))
419	>	end if
420
421	<	return
422	<	end subroutine setSimVariables
421	>	groupMaxCutoffCol = 0.0_dp
422	>	gtypeMaxCutoffCol = 0.0_dp
423
424	+	#endif
425	+	groupMaxCutoffRow = 0.0_dp
426	+	gtypeMaxCutoffRow = 0.0_dp
427	+
428	+
429	+	!! first we do a single loop over the cutoff groups to find the
430	+	!! largest cutoff for any atypes present in this group.  We also
431	+	!! create gtypes at this point.
432	+
433	+	tol = 1.0d-6
434	+	nGroupTypesRow = 0
435	+	nGroupTypesCol = 0
436	+	do i = istart, iend
437	+	n_in_i = groupStartRow(i+1) - groupStartRow(i)
438	+	groupMaxCutoffRow(i) = 0.0_dp
439	+	do ia = groupStartRow(i), groupStartRow(i+1)-1
440	+	atom1 = groupListRow(ia)
441	+	#ifdef IS_MPI
442	+	me_i = atid_row(atom1)
443	+	#else
444	+	me_i = atid(atom1)
445	+	#endif
446	+	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
447	+	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
448	+	endif
449	+	enddo
450	+	if (nGroupTypesRow.eq.0) then
451	+	nGroupTypesRow = nGroupTypesRow + 1
452	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
453	+	groupToGtypeRow(i) = nGroupTypesRow
454	+	else
455	+	GtypeFound = .false.
456	+	do g = 1, nGroupTypesRow
457	+	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
458	+	groupToGtypeRow(i) = g
459	+	GtypeFound = .true.
460	+	endif
461	+	enddo
462	+	if (.not.GtypeFound) then
463	+	nGroupTypesRow = nGroupTypesRow + 1
464	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
465	+	groupToGtypeRow(i) = nGroupTypesRow
466	+	endif
467	+	endif
468	+	enddo
469	+
470	+	#ifdef IS_MPI
471	+	do j = jstart, jend
472	+	n_in_j = groupStartCol(j+1) - groupStartCol(j)
473	+	groupMaxCutoffCol(j) = 0.0_dp
474	+	do ja = groupStartCol(j), groupStartCol(j+1)-1
475	+	atom1 = groupListCol(ja)
476	+
477	+	me_j = atid_col(atom1)
478	+
479	+	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
480	+	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
481	+	endif
482	+	enddo
483	+
484	+	if (nGroupTypesCol.eq.0) then
485	+	nGroupTypesCol = nGroupTypesCol + 1
486	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
487	+	groupToGtypeCol(j) = nGroupTypesCol
488	+	else
489	+	GtypeFound = .false.
490	+	do g = 1, nGroupTypesCol
491	+	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
492	+	groupToGtypeCol(j) = g
493	+	GtypeFound = .true.
494	+	endif
495	+	enddo
496	+	if (.not.GtypeFound) then
497	+	nGroupTypesCol = nGroupTypesCol + 1
498	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
499	+	groupToGtypeCol(j) = nGroupTypesCol
500	+	endif
501	+	endif
502	+	enddo
503	+
504	+	#else
505	+	! Set pointers to information we just found
506	+	nGroupTypesCol = nGroupTypesRow
507	+	groupToGtypeCol => groupToGtypeRow
508	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
509	+	groupMaxCutoffCol => groupMaxCutoffRow
510	+	#endif
511	+
512	+	!! allocate the gtypeCutoffMap here.
513	+	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
514	+	!! then we do a double loop over all the group TYPES to find the cutoff
515	+	!! map between groups of two types
516	+	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
517	+
518	+	do i = 1, nGroupTypesRow
519	+	do j = 1, nGroupTypesCol
520	+
521	+	select case(cutoffPolicy)
522	+	case(TRADITIONAL_CUTOFF_POLICY)
523	+	thisRcut = tradRcut
524	+	case(MIX_CUTOFF_POLICY)
525	+	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
526	+	case(MAX_CUTOFF_POLICY)
527	+	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
528	+	case default
529	+	call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
530	+	return
531	+	end select
532	+	gtypeCutoffMap(i,j)%rcut = thisRcut
533	+
534	+	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
535	+
536	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
537	+
538	+	if (.not.haveSkinThickness) then
539	+	skinThickness = 1.0_dp
540	+	endif
541	+
542	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
543	+
544	+	! sanity check
545	+
546	+	if (haveDefaultCutoffs) then
547	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
548	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
549	+	endif
550	+	endif
551	+	enddo
552	+	enddo
553	+
554	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
555	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
556	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
557	+	#ifdef IS_MPI
558	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
559	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
560	+	#endif
561	+	groupMaxCutoffCol => null()
562	+	gtypeMaxCutoffCol => null()
563	+
564	+	haveGtypeCutoffMap = .true.
565	+	end subroutine createGtypeCutoffMap
566	+
567	+	subroutine setCutoffs(defRcut, defRsw)
568	+
569	+	real(kind=dp),intent(in) :: defRcut, defRsw
570	+	character(len = statusMsgSize) :: errMsg
571	+	integer :: localError
572	+
573	+	defaultRcut = defRcut
574	+	defaultRsw = defRsw
575	+
576	+	defaultDoShift = .false.
577	+	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
578	+
579	+	write(errMsg, *) &
580	+	'cutoffRadius and switchingRadius are set to the same', newline &
581	+	// tab, 'value.  OOPSE will use shifted ', newline &
582	+	// tab, 'potentials instead of switching functions.'
583	+
584	+	call handleInfo("setCutoffs", errMsg)
585	+
586	+	defaultDoShift = .true.
587	+
588	+	endif
589	+
590	+	localError = 0
591	+	call setLJDefaultCutoff( defaultRcut, defaultDoShift )
592	+	call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
593	+	call setCutoffEAM( defaultRcut )
594	+	call setCutoffSC( defaultRcut )
595	+	call set_switch(defaultRsw, defaultRcut)
596	+	call setHmatDangerousRcutValue(defaultRcut)
597	+	call setupSplineSqrt(defaultRcut)
598	+
599	+	haveDefaultCutoffs = .true.
600	+	haveGtypeCutoffMap = .false.
601	+
602	+	end subroutine setCutoffs
603	+
604	+	subroutine cWasLame()
605	+
606	+	VisitCutoffsAfterComputing = .true.
607	+	return
608	+
609	+	end subroutine cWasLame
610	+
611	+	subroutine setCutoffPolicy(cutPolicy)
612	+
613	+	integer, intent(in) :: cutPolicy
614	+
615	+	cutoffPolicy = cutPolicy
616	+	haveCutoffPolicy = .true.
617	+	haveGtypeCutoffMap = .false.
618	+
619	+	end subroutine setCutoffPolicy
620	+
621	+	subroutine setElectrostaticMethod( thisESM )
622	+
623	+	integer, intent(in) :: thisESM
624	+
625	+	electrostaticSummationMethod = thisESM
626	+	haveElectrostaticSummationMethod = .true.
627	+
628	+	end subroutine setElectrostaticMethod
629	+
630	+	subroutine setSkinThickness( thisSkin )
631	+
632	+	real(kind=dp), intent(in) :: thisSkin
633	+
634	+	skinThickness = thisSkin
635	+	haveSkinThickness = .true.
636	+	haveGtypeCutoffMap = .false.
637	+
638	+	end subroutine setSkinThickness
639	+
640	+	subroutine setSimVariables()
641	+	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
642	+	SIM_uses_EAM = SimUsesEAM()
643	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
644	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
645	+	SIM_uses_PBC = SimUsesPBC()
646	+	SIM_uses_SC = SimUsesSC()
647	+
648	+	haveSIMvariables = .true.
649	+
650	+	return
651	+	end subroutine setSimVariables
652	+
653		subroutine doReadyCheck(error)
654		integer, intent(out) :: error
380	–
655		integer :: myStatus
656
657		error = 0
658
659	<	if (.not. haveInteractionMap) then
659	>	if (.not. haveInteractionHash) then
660	>	call createInteractionHash()
661	>	endif
662
663	<	myStatus = 0
663	>	if (.not. haveGtypeCutoffMap) then
664	>	call createGtypeCutoffMap()
665	>	endif
666
667	<	call createInteractionMap(myStatus)
667	>	if (VisitCutoffsAfterComputing) then
668	>	call set_switch(largestRcut, largestRcut)
669	>	call setHmatDangerousRcutValue(largestRcut)
670	>	call setLJsplineRmax(largestRcut)
671	>	call setCutoffEAM(largestRcut)
672	>	call setCutoffSC(largestRcut)
673	>	VisitCutoffsAfterComputing = .false.
674	>	endif
675
676	<	if (myStatus .ne. 0) then
677	<	write(default_error, *) 'createInteractionMap failed in doForces!'
393	<	error = -1
394	<	return
395	<	endif
676	>	if (.not. haveSplineSqrt) then
677	>	call setupSplineSqrt(largestRcut)
678		endif
679
680		if (.not. haveSIMvariables) then
681		call setSimVariables()
682		endif
683
402	–	if (.not. haveRlist) then
403	–	write(default_error, *) 'rList has not been set in doForces!'
404	–	error = -1
405	–	return
406	–	endif
407	–
684		if (.not. haveNeighborList) then
685		write(default_error, *) 'neighbor list has not been initialized in doForces!'
686		error = -1
687		return
688		end if
689	<
689	>
690		if (.not. haveSaneForceField) then
691		write(default_error, *) 'Force Field is not sane in doForces!'
692		error = -1
693		return
694		end if
695	<
695	>
696		#ifdef IS_MPI
697		if (.not. isMPISimSet()) then
698		write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
#	Line 428 | Line 704 | contains
704		end subroutine doReadyCheck
705
706
707	<	subroutine init_FF(use_RF_c, thisStat)
707	>	subroutine init_FF(thisStat)
708
433	–	logical, intent(in) :: use_RF_c
434	–
709		integer, intent(out) :: thisStat
710		integer :: my_status, nMatches
711		integer, pointer :: MatchList(:) => null()
438	–	real(kind=dp) :: rcut, rrf, rt, dielect
712
713		!! assume things are copacetic, unless they aren't
714		thisStat = 0
715
443	–	!! Fortran's version of a cast:
444	–	FF_uses_RF = use_RF_c
445	–
716		!! init_FF is called *after* all of the atom types have been
717		!! defined in atype_module using the new_atype subroutine.
718		!!
#	Line 450 | Line 720 | contains
720		!! interactions are used by the force field.
721
722		FF_uses_DirectionalAtoms = .false.
453	–	FF_uses_LennardJones = .false.
454	–	FF_uses_Electrostatics = .false.
455	–	FF_uses_Charges = .false.
723		FF_uses_Dipoles = .false.
457	–	FF_uses_Sticky = .false.
458	–	FF_uses_StickyPower = .false.
724		FF_uses_GayBerne = .false.
725		FF_uses_EAM = .false.
726	<	FF_uses_Shapes = .false.
462	<	FF_uses_FLARB = .false.
726	>	FF_uses_SC = .false.
727
728		call getMatchingElementList(atypes, "is_Directional", .true., &
729		nMatches, MatchList)
730		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
731
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	–
732		call getMatchingElementList(atypes, "is_Dipole", .true., &
733		nMatches, MatchList)
734	<	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
734	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
735
736		call getMatchingElementList(atypes, "is_GayBerne", .true., &
737		nMatches, MatchList)
738	<	if (nMatches .gt. 0) then
518	<	FF_uses_GayBerne = .true.
519	<	FF_uses_DirectionalAtoms = .true.
520	<	endif
738	>	if (nMatches .gt. 0) FF_uses_GayBerne = .true.
739
740		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
741		if (nMatches .gt. 0) FF_uses_EAM = .true.
742
743	<	call getMatchingElementList(atypes, "is_Shape", .true., &
744	<	nMatches, MatchList)
527	<	if (nMatches .gt. 0) then
528	<	FF_uses_Shapes = .true.
529	<	FF_uses_DirectionalAtoms = .true.
530	<	endif
743	>	call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
744	>	if (nMatches .gt. 0) FF_uses_SC = .true.
745
532	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
533	–	nMatches, MatchList)
534	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
746
536	–	!! Assume sanity (for the sake of argument)
747		haveSaneForceField = .true.
538	–
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
748
749		if (FF_uses_EAM) then
750		call init_EAM_FF(my_status)
#	Line 572 | Line 756 | contains
756		end if
757		endif
758
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	–
759		if (.not. haveNeighborList) then
760		!! Create neighbor lists
761		call expandNeighborList(nLocal, my_status)
#	Line 617 | Line 789 | contains
789
790		!! Stress Tensor
791		real( kind = dp), dimension(9) :: tau
792	<	real ( kind = dp ) :: pot
792	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
793		logical ( kind = 2) :: do_pot_c, do_stress_c
794		logical :: do_pot
795		logical :: do_stress
796		logical :: in_switching_region
797		#ifdef IS_MPI
798	<	real( kind = DP ) :: pot_local
798	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
799		integer :: nAtomsInRow
800		integer :: nAtomsInCol
801		integer :: nprocs
#	Line 638 | Line 810 | contains
810		integer :: nlist
811		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
812		real( kind = DP ) :: sw, dswdr, swderiv, mf
813	+	real( kind = DP ) :: rVal
814		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
815		real(kind=dp) :: rfpot, mu_i, virial
816	+	real(kind=dp):: rCut
817		integer :: me_i, me_j, n_in_i, n_in_j
818		logical :: is_dp_i
819		integer :: neighborListSize
#	Line 647 | Line 821 | contains
821		integer :: localError
822		integer :: propPack_i, propPack_j
823		integer :: loopStart, loopEnd, loop
824	<	integer :: iMap
825	<	real(kind=dp) :: listSkin = 1.0
824	>	integer :: iHash
825	>	integer :: i1
826	>
827
828		!! initialize local variables
829
#	Line 712 | Line 887 | contains
887		! (but only on the first time through):
888		if (loop .eq. loopStart) then
889		#ifdef IS_MPI
890	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
890	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
891		update_nlist)
892		#else
893	<	call checkNeighborList(nGroups, q_group, listSkin, &
893	>	call checkNeighborList(nGroups, q_group, skinThickness, &
894		update_nlist)
895		#endif
896		endif
#	Line 739 | Line 914 | contains
914		#endif
915		outer: do i = istart, iend
916
742	–	#ifdef IS_MPI
743	–	me_i = atid_row(i)
744	–	#else
745	–	me_i = atid(i)
746	–	#endif
747	–
917		if (update_nlist) point(i) = nlist + 1
918
919		n_in_i = groupStartRow(i+1) - groupStartRow(i)
#	Line 779 | Line 948 | contains
948		me_j = atid(j)
949		call get_interatomic_vector(q_group(:,i), &
950		q_group(:,j), d_grp, rgrpsq)
951	<	#endif
951	>	#endif
952
953	<	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
953	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
954		if (update_nlist) then
955		nlist = nlist + 1
956
#	Line 801 | Line 970 | contains
970
971		list(nlist) = j
972		endif
973	<
974	<	if (loop .eq. PAIR_LOOP) then
806	<	vij = 0.0d0
807	<	fij(1:3) = 0.0d0
808	<	endif
973	>
974	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
975
976	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
977	<	in_switching_region)
978	<
979	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
980	<
981	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
982	<
983	<	atom1 = groupListRow(ia)
984	<
985	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
986	<
987	<	atom2 = groupListCol(jb)
988	<
989	<	if (skipThisPair(atom1, atom2)) cycle inner
990	<
991	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
992	<	d_atm(1:3) = d_grp(1:3)
993	<	ratmsq = rgrpsq
994	<	else
976	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
977	>	if (loop .eq. PAIR_LOOP) then
978	>	vij = 0.0d0
979	>	fij(1) = 0.0_dp
980	>	fij(2) = 0.0_dp
981	>	fij(3) = 0.0_dp
982	>	endif
983	>
984	>	call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
985	>
986	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
987	>
988	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
989	>
990	>	atom1 = groupListRow(ia)
991	>
992	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
993	>
994	>	atom2 = groupListCol(jb)
995	>
996	>	if (skipThisPair(atom1, atom2))  cycle inner
997	>
998	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
999	>	d_atm(1) = d_grp(1)
1000	>	d_atm(2) = d_grp(2)
1001	>	d_atm(3) = d_grp(3)
1002	>	ratmsq = rgrpsq
1003	>	else
1004		#ifdef IS_MPI
1005	<	call get_interatomic_vector(q_Row(:,atom1), &
1006	<	q_Col(:,atom2), d_atm, ratmsq)
1005	>	call get_interatomic_vector(q_Row(:,atom1), &
1006	>	q_Col(:,atom2), d_atm, ratmsq)
1007		#else
1008	<	call get_interatomic_vector(q(:,atom1), &
1009	<	q(:,atom2), d_atm, ratmsq)
1008	>	call get_interatomic_vector(q(:,atom1), &
1009	>	q(:,atom2), d_atm, ratmsq)
1010		#endif
1011	<	endif
1012	<
1013	<	if (loop .eq. PREPAIR_LOOP) then
1011	>	endif
1012	>
1013	>	if (loop .eq. PREPAIR_LOOP) then
1014		#ifdef IS_MPI
1015	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1016	<	rgrpsq, d_grp, do_pot, do_stress, &
1017	<	eFrame, A, f, t, pot_local)
1015	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1016	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1017	>	eFrame, A, f, t, pot_local)
1018		#else
1019	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1020	<	rgrpsq, d_grp, do_pot, do_stress, &
1021	<	eFrame, A, f, t, pot)
1019	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1020	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1021	>	eFrame, A, f, t, pot)
1022		#endif
1023	<	else
1023	>	else
1024		#ifdef IS_MPI
1025	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1026	<	do_pot, &
1027	<	eFrame, A, f, t, pot_local, vpair, fpair)
1025	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1026	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1027	>	fpair, d_grp, rgrp, rCut)
1028		#else
1029	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1030	<	do_pot,  &
1031	<	eFrame, A, f, t, pot, vpair, fpair)
1029	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1030	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1031	>	d_grp, rgrp, rCut)
1032		#endif
1033	+	vij = vij + vpair
1034	+	fij(1) = fij(1) + fpair(1)
1035	+	fij(2) = fij(2) + fpair(2)
1036	+	fij(3) = fij(3) + fpair(3)
1037	+	endif
1038	+	enddo inner
1039	+	enddo
1040
1041	<	vij = vij + vpair
1042	<	fij(1:3) = fij(1:3) + fpair(1:3)
1043	<	endif
1044	<	enddo inner
1045	<	enddo
1046	<
1047	<	if (loop .eq. PAIR_LOOP) then
1048	<	if (in_switching_region) then
1049	<	swderiv = vij*dswdr/rgrp
1050	<	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)
1041	>	if (loop .eq. PAIR_LOOP) then
1042	>	if (in_switching_region) then
1043	>	swderiv = vij*dswdr/rgrp
1044	>	fij(1) = fij(1) + swderiv*d_grp(1)
1045	>	fij(2) = fij(2) + swderiv*d_grp(2)
1046	>	fij(3) = fij(3) + swderiv*d_grp(3)
1047	>
1048	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1049	>	atom1=groupListRow(ia)
1050	>	mf = mfactRow(atom1)
1051		#ifdef IS_MPI
1052	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1053	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1054	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1052	>	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1053	>	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1054	>	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1055		#else
1056	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1057	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1058	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1056	>	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1057	>	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1058	>	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1059		#endif
1060	<	enddo
1061	<
1062	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1063	<	atom2=groupListCol(jb)
1064	<	mf = mfactCol(atom2)
1060	>	enddo
1061	>
1062	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1063	>	atom2=groupListCol(jb)
1064	>	mf = mfactCol(atom2)
1065		#ifdef IS_MPI
1066	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1067	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1068	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1066	>	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1067	>	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1068	>	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1069		#else
1070	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1071	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1072	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1070	>	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1071	>	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1072	>	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1073		#endif
1074	<	enddo
1075	<	endif
1074	>	enddo
1075	>	endif
1076
1077	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1077	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1078	>	endif
1079		endif
1080	<	end if
1080	>	endif
1081		enddo
1082	+
1083		enddo outer
1084
1085		if (update_nlist) then
#	Line 961 | Line 1139 | contains
1139
1140		if (do_pot) then
1141		! scatter/gather pot_row into the members of my column
1142	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1143	<
1142	>	do i = 1,LR_POT_TYPES
1143	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1144	>	end do
1145		! scatter/gather pot_local into all other procs
1146		! add resultant to get total pot
1147		do i = 1, nlocal
1148	<	pot_local = pot_local + pot_Temp(i)
1148	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1149	>	+ pot_Temp(1:LR_POT_TYPES,i)
1150		enddo
1151
1152		pot_Temp = 0.0_DP
1153	<
1154	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1153	>	do i = 1,LR_POT_TYPES
1154	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1155	>	end do
1156		do i = 1, nlocal
1157	<	pot_local = pot_local + pot_Temp(i)
1157	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1158	>	+ pot_Temp(1:LR_POT_TYPES,i)
1159		enddo
1160
1161		endif
1162		#endif
1163
1164	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1164	>	if (SIM_requires_postpair_calc) then
1165	>	do i = 1, nlocal
1166	>
1167	>	! we loop only over the local atoms, so we don't need row and column
1168	>	! lookups for the types
1169	>
1170	>	me_i = atid(i)
1171	>
1172	>	! is the atom electrostatic?  See if it would have an
1173	>	! electrostatic interaction with itself
1174	>	iHash = InteractionHash(me_i,me_i)
1175
1176	<	if (FF_uses_RF .and. SIM_uses_RF) then
985	<
1176	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1177		#ifdef IS_MPI
1178	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1179	<	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)
1178	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1179	>	t, do_pot)
1180		#else
1181	<	me_i = atid(i)
1181	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1182	>	t, do_pot)
1183		#endif
1184	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1185	<
1186	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1187	<
1188	<	mu_i = getDipoleMoment(me_i)
1189	<
1190	<	!! The reaction field needs to include a self contribution
1191	<	!! to the field:
1192	<	call accumulate_self_rf(i, mu_i, eFrame)
1193	<	!! Get the reaction field contribution to the
1194	<	!! potential and torques:
1195	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1196	<	#ifdef IS_MPI
1197	<	pot_local = pot_local + rfpot
1184	>	endif
1185	>
1186	>
1187	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1188	>
1189	>	! loop over the excludes to accumulate RF stuff we've
1190	>	! left out of the normal pair loop
1191	>
1192	>	do i1 = 1, nSkipsForAtom(i)
1193	>	j = skipsForAtom(i, i1)
1194	>
1195	>	! prevent overcounting of the skips
1196	>	if (i.lt.j) then
1197	>	call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1198	>	rVal = dsqrt(ratmsq)
1199	>	call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1200	>	#ifdef IS_MPI
1201	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1202	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1203		#else
1204	<	pot = pot + rfpot
1205	<
1204	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1205	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1206		#endif
1207	<	endif
1208	<	enddo
1209	<	endif
1207	>	endif
1208	>	enddo
1209	>	endif
1210	>	enddo
1211		endif
1212	<
1025	<
1212	>
1213		#ifdef IS_MPI
1214	<
1214	>
1215		if (do_pot) then
1216	<	pot = pot + pot_local
1217	<	!! we assume the c code will do the allreduce to get the total potential
1031	<	!! we could do it right here if we needed to...
1216	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1217	>	mpi_comm_world,mpi_err)
1218		endif
1219	<
1219	>
1220		if (do_stress) then
1221		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1222		mpi_comm_world,mpi_err)
1223		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1224		mpi_comm_world,mpi_err)
1225		endif
1226	<
1226	>
1227		#else
1228	<
1228	>
1229		if (do_stress) then
1230		tau = tau_Temp
1231		virial = virial_Temp
1232		endif
1233	<
1233	>
1234		#endif
1235	<
1235	>
1236		end subroutine do_force_loop
1237
1238		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1239	<	eFrame, A, f, t, pot, vpair, fpair)
1239	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1240
1241	<	real( kind = dp ) :: pot, vpair, sw
1241	>	real( kind = dp ) :: vpair, sw
1242	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1243		real( kind = dp ), dimension(3) :: fpair
1244		real( kind = dp ), dimension(nLocal)   :: mfact
1245		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1063 | Line 1250 | contains
1250		logical, intent(inout) :: do_pot
1251		integer, intent(in) :: i, j
1252		real ( kind = dp ), intent(inout) :: rijsq
1253	<	real ( kind = dp )                :: r
1253	>	real ( kind = dp ), intent(inout) :: r_grp
1254		real ( kind = dp ), intent(inout) :: d(3)
1255	<	real ( kind = dp ) :: ebalance
1255	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1256	>	real ( kind = dp ), intent(inout) :: rCut
1257	>	real ( kind = dp ) :: r
1258	>	real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1259		integer :: me_i, me_j
1260	+	integer :: k
1261
1262	<	integer :: iMap
1262	>	integer :: iHash
1263
1264	<	r = sqrt(rijsq)
1264	>	if (useSpline) then
1265	>	call lookupUniformSpline(splineSqrt, rijsq, r)
1266	>	else
1267	>	r = sqrt(rijsq)
1268	>	endif
1269	>
1270		vpair = 0.0d0
1271		fpair(1:3) = 0.0d0
1272
#	Line 1082 | Line 1278 | contains
1278		me_j = atid(j)
1279		#endif
1280
1281	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1282	<
1283	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1284	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1281	>	iHash = InteractionHash(me_i, me_j)
1282	>
1283	>	if ( iand(iHash, LJ_PAIR).ne.0 ) then
1284	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1285	>	pot(VDW_POT), f, do_pot)
1286		endif
1287	<
1288	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1289	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1290	<	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	<
1287	>
1288	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1289	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1290	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1291		endif
1292	<
1293	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1292	>
1293	>	if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1294		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1295	<	pot, A, f, t, do_pot)
1295	>	pot(HB_POT), A, f, t, do_pot)
1296		endif
1297	<
1298	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1297	>
1298	>	if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1299		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1300	<	pot, A, f, t, do_pot)
1300	>	pot(HB_POT), A, f, t, do_pot)
1301		endif
1302	<
1303	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1302	>
1303	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1304		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1305	<	pot, A, f, t, do_pot)
1305	>	pot(VDW_POT), A, f, t, do_pot)
1306		endif
1307
1308	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1309	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1310	<	!           pot, A, f, t, do_pot)
1308	>	if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1309	>	call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1310	>	pot(VDW_POT), A, f, t, do_pot)
1311		endif
1312	<
1313	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1314	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1315	<	do_pot)
1312	>
1313	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1314	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1315	>	pot(METALLIC_POT), f, do_pot)
1316		endif
1317	<
1318	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1317	>
1318	>	if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1319		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1320	<	pot, A, f, t, do_pot)
1320	>	pot(VDW_POT), A, f, t, do_pot)
1321		endif
1322	<
1323	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1322	>
1323	>	if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1324		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1325	<	pot, A, f, t, do_pot)
1325	>	pot(VDW_POT), A, f, t, do_pot)
1326	>	endif
1327	>
1328	>	if ( iand(iHash, SC_PAIR).ne.0 ) then
1329	>	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1330	>	pot(METALLIC_POT), f, do_pot)
1331		endif
1332	<
1332	>
1333		end subroutine do_pair
1334
1335	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1335	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1336		do_pot, do_stress, eFrame, A, f, t, pot)
1337
1338	<	real( kind = dp ) :: pot, sw
1338	>	real( kind = dp ) :: sw
1339	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1340		real( kind = dp ), dimension(9,nLocal) :: eFrame
1341		real (kind=dp), dimension(9,nLocal) :: A
1342		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1149 | Line 1344 | contains
1344
1345		logical, intent(inout) :: do_pot, do_stress
1346		integer, intent(in) :: i, j
1347	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1347	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1348		real ( kind = dp )                :: r, rc
1349		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1350
1351	<	integer :: me_i, me_j, iMap
1351	>	integer :: me_i, me_j, iHash
1352
1353	+	if (useSpline) then
1354	+	call lookupUniformSpline(splineSqrt, rijsq, r)
1355	+	else
1356	+	r = sqrt(rijsq)
1357	+	endif
1358	+
1359		#ifdef IS_MPI
1360		me_i = atid_row(i)
1361		me_j = atid_col(j)
#	Line 1163 | Line 1364 | contains
1364		me_j = atid(j)
1365		#endif
1366
1367	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1367	>	iHash = InteractionHash(me_i, me_j)
1368
1369	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1370	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1369	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1370	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1371		endif
1372	+
1373	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1374	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1375	+	endif
1376
1377		end subroutine do_prepair
1378
1379
1380		subroutine do_preforce(nlocal,pot)
1381		integer :: nlocal
1382	<	real( kind = dp ) :: pot
1382	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1383
1384		if (FF_uses_EAM .and. SIM_uses_EAM) then
1385	<	call calc_EAM_preforce_Frho(nlocal,pot)
1385	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1386		endif
1387	<
1388	<
1387	>	if (FF_uses_SC .and. SIM_uses_SC) then
1388	>	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1389	>	endif
1390		end subroutine do_preforce
1391
1392
#	Line 1192 | Line 1398 | contains
1398		real( kind = dp ) :: d(3), scaled(3)
1399		integer i
1400
1401	<	d(1:3) = q_j(1:3) - q_i(1:3)
1401	>	d(1) = q_j(1) - q_i(1)
1402	>	d(2) = q_j(2) - q_i(2)
1403	>	d(3) = q_j(3) - q_i(3)
1404
1405		! Wrap back into periodic box if necessary
1406		if ( SIM_uses_PBC ) then
1407
1408		if( .not.boxIsOrthorhombic ) then
1409		! calc the scaled coordinates.
1410	+	! scaled = matmul(HmatInv, d)
1411
1412	<	scaled = matmul(HmatInv, d)
1413	<
1412	>	scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1413	>	scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1414	>	scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1415	>
1416		! wrap the scaled coordinates
1417
1418	<	scaled = scaled  - anint(scaled)
1418	>	scaled(1) = scaled(1) - dnint(scaled(1))
1419	>	scaled(2) = scaled(2) - dnint(scaled(2))
1420	>	scaled(3) = scaled(3) - dnint(scaled(3))
1421
1209	–
1422		! calc the wrapped real coordinates from the wrapped scaled
1423		! coordinates
1424	+	! d = matmul(Hmat,scaled)
1425	+	d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1426	+	d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1427	+	d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1428
1213	–	d = matmul(Hmat,scaled)
1214	–
1429		else
1430		! calc the scaled coordinates.
1431
1432	<	do i = 1, 3
1433	<	scaled(i) = d(i) * HmatInv(i,i)
1434	<
1435	<	! wrap the scaled coordinates
1436	<
1437	<	scaled(i) = scaled(i) - anint(scaled(i))
1432	>	scaled(1) = d(1) * HmatInv(1,1)
1433	>	scaled(2) = d(2) * HmatInv(2,2)
1434	>	scaled(3) = d(3) * HmatInv(3,3)
1435	>
1436	>	! wrap the scaled coordinates
1437	>
1438	>	scaled(1) = scaled(1) - dnint(scaled(1))
1439	>	scaled(2) = scaled(2) - dnint(scaled(2))
1440	>	scaled(3) = scaled(3) - dnint(scaled(3))
1441
1442	<	! calc the wrapped real coordinates from the wrapped scaled
1443	<	! coordinates
1442	>	! calc the wrapped real coordinates from the wrapped scaled
1443	>	! coordinates
1444
1445	<	d(i) = scaled(i)*Hmat(i,i)
1446	<	enddo
1445	>	d(1) = scaled(1)*Hmat(1,1)
1446	>	d(2) = scaled(2)*Hmat(2,2)
1447	>	d(3) = scaled(3)*Hmat(3,3)
1448	>
1449		endif
1450
1451		endif
1452
1453	<	r_sq = dot_product(d,d)
1453	>	r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1454
1455		end subroutine get_interatomic_vector
1456
#	Line 1263 | Line 1482 | contains
1482		pot_Col = 0.0_dp
1483		pot_Temp = 0.0_dp
1484
1266	–	rf_Row = 0.0_dp
1267	–	rf_Col = 0.0_dp
1268	–	rf_Temp = 0.0_dp
1269	–
1485		#endif
1486
1487		if (FF_uses_EAM .and. SIM_uses_EAM) then
1488		call clean_EAM()
1489		endif
1490
1276	–	rf = 0.0_dp
1491		tau_Temp = 0.0_dp
1492		virial_Temp = 0.0_dp
1493		end subroutine zero_work_arrays
#	Line 1362 | Line 1576 | contains
1576
1577		function FF_UsesDirectionalAtoms() result(doesit)
1578		logical :: doesit
1579	<	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
1579	>	doesit = FF_uses_DirectionalAtoms
1580		end function FF_UsesDirectionalAtoms
1581
1582		function FF_RequiresPrepairCalc() result(doesit)
1583		logical :: doesit
1584	<	doesit = FF_uses_EAM
1584	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1585	>	.or. FF_uses_MEAM
1586		end function FF_RequiresPrepairCalc
1587
1375	–	function FF_RequiresPostpairCalc() result(doesit)
1376	–	logical :: doesit
1377	–	doesit = FF_uses_RF
1378	–	end function FF_RequiresPostpairCalc
1379	–
1588		#ifdef PROFILE
1589		function getforcetime() result(totalforcetime)
1590		real(kind=dp) :: totalforcetime
#	Line 1409 | Line 1617 | end module doForces
1617
1618		end subroutine add_stress_tensor
1619
1620	+	subroutine setupSplineSqrt(rmax)
1621	+	real(kind=dp), intent(in) :: rmax
1622	+	real(kind=dp), dimension(np) :: xvals, yvals
1623	+	real(kind=dp) :: r2_1, r2_n, dx, r2
1624	+	integer :: i
1625	+
1626	+	r2_1 = 0.5d0
1627	+	r2_n = rmax*rmax
1628	+
1629	+	dx = (r2_n-r2_1) / dble(np-1)
1630	+
1631	+	do i = 1, np
1632	+	r2 = r2_1 + dble(i-1)*dx
1633	+	xvals(i) = r2
1634	+	yvals(i) = dsqrt(r2)
1635	+	enddo
1636	+
1637	+	call newSpline(splineSqrt, xvals, yvals, .true.)
1638	+
1639	+	haveSplineSqrt = .true.
1640	+	return
1641	+	end subroutine setupSplineSqrt
1642	+
1643	+	subroutine deleteSplineSqrt()
1644	+	call deleteSpline(splineSqrt)
1645	+	haveSplineSqrt = .false.
1646	+	return
1647	+	end subroutine deleteSplineSqrt
1648	+
1649		end module doForces

Diff Legend

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