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 2262 by chuckv, Sun Jul 3 20:53:43 2005 UTC vs.
Revision 2787 by gezelter, Mon Jun 5 18:24:45 2006 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.23 2005-07-03 20:53:43 chuckv Exp $, $Date: 2005-07-03 20:53:43 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $
48	>	!! @version $Id: doForces.F90,v 1.83 2006-06-05 18:24:45 gezelter Exp $, $Date: 2006-06-05 18:24:45 $, $Name: not supported by cvs2svn $, $Revision: 1.83 $
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
93		logical, save :: FF_uses_DirectionalAtoms
88	–	logical, save :: FF_uses_LennardJones
89	–	logical, save :: FF_uses_Electrostatics
90	–	logical, save :: FF_uses_Charges
94		logical, save :: FF_uses_Dipoles
92	–	logical, save :: FF_uses_Quadrupoles
93	–	logical, save :: FF_uses_Sticky
94	–	logical, save :: FF_uses_StickyPower
95		logical, save :: FF_uses_GayBerne
96		logical, save :: FF_uses_EAM
97	<	logical, save :: FF_uses_Shapes
98	<	logical, save :: FF_uses_FLARB
99	<	logical, save :: FF_uses_RF
97	>	logical, save :: FF_uses_SC
98	>	logical, save :: FF_uses_MEAM
99	>
100
101		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
102		logical, save :: SIM_uses_EAM
103	<	logical, save :: SIM_uses_Shapes
104	<	logical, save :: SIM_uses_FLARB
113	<	logical, save :: SIM_uses_RF
103	>	logical, save :: SIM_uses_SC
104	>	logical, save :: SIM_uses_MEAM
105		logical, save :: SIM_requires_postpair_calc
106		logical, save :: SIM_requires_prepair_calc
107		logical, save :: SIM_uses_PBC
117	–	logical, save :: SIM_uses_molecular_cutoffs
108
109	<	!!!GO AWAY---------
110	<	!!!!!real(kind=dp), save :: rlist, rlistsq
109	>	integer, save :: electrostaticSummationMethod
110	>	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
111
112	+	real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
113	+	real(kind=dp), save :: skinThickness
114	+	logical, save :: defaultDoShift
115	+
116		public :: init_FF
117	+	public :: setCutoffs
118	+	public :: cWasLame
119	+	public :: setElectrostaticMethod
120	+	public :: setCutoffPolicy
121	+	public :: setSkinThickness
122		public :: do_force_loop
124	–	!  public :: setRlistDF
125	–	!public :: addInteraction
126	–	!public :: setInteractionHash
127	–	!public :: getInteractionHash
128	–	public :: createInteractionMap
129	–	public :: createRcuts
123
124		#ifdef PROFILE
125		public :: getforcetime
#	Line 134 | Line 127 | module doForces
127		real :: forceTimeInitial, forceTimeFinal
128		integer :: nLoops
129		#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
130
131	<	type(Interaction), dimension(:,:),allocatable :: InteractionMap
132	<
131	>	!! Variables for cutoff mapping and interaction mapping
132	>	! Bit hash to determine pair-pair interactions.
133	>	integer, dimension(:,:), allocatable :: InteractionHash
134	>	real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
135	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
136	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
137
138	<
138	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
139	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
140	>
141	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
142	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
143	>	type ::gtypeCutoffs
144	>	real(kind=dp) :: rcut
145	>	real(kind=dp) :: rcutsq
146	>	real(kind=dp) :: rlistsq
147	>	end type gtypeCutoffs
148	>	type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
149	>
150	>
151		contains
152
153	<
151	<	subroutine createInteractionMap(status)
153	>	subroutine createInteractionHash()
154		integer :: nAtypes
153	–	integer :: status
155		integer :: i
156		integer :: j
157	<	integer :: ihash
158	<	real(kind=dp) :: myRcut
158	<	! Test Types
157	>	integer :: iHash
158	>	!! Test Types
159		logical :: i_is_LJ
160		logical :: i_is_Elect
161		logical :: i_is_Sticky
#	Line 163 | Line 163 | contains
163		logical :: i_is_GB
164		logical :: i_is_EAM
165		logical :: i_is_Shape
166	+	logical :: i_is_SC
167	+	logical :: i_is_MEAM
168		logical :: j_is_LJ
169		logical :: j_is_Elect
170		logical :: j_is_Sticky
#	Line 170 | Line 172 | contains
172		logical :: j_is_GB
173		logical :: j_is_EAM
174		logical :: j_is_Shape
175	<
176	<
175	>	logical :: j_is_SC
176	>	logical :: j_is_MEAM
177	>	real(kind=dp) :: myRcut
178	>
179		if (.not. associated(atypes)) then
180	<	call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
177	<	status = -1
180	>	call handleError("doForces", "atypes was not present before call of createInteractionHash!")
181		return
182		endif
183
184		nAtypes = getSize(atypes)
185
186		if (nAtypes == 0) then
187	<	status = -1
187	>	call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
188		return
189		end if
190
191	<	if (.not. allocated(InteractionMap)) then
192	<	allocate(InteractionMap(nAtypes,nAtypes))
191	>	if (.not. allocated(InteractionHash)) then
192	>	allocate(InteractionHash(nAtypes,nAtypes))
193	>	else
194	>	deallocate(InteractionHash)
195	>	allocate(InteractionHash(nAtypes,nAtypes))
196		endif
197	+
198	+	if (.not. allocated(atypeMaxCutoff)) then
199	+	allocate(atypeMaxCutoff(nAtypes))
200	+	else
201	+	deallocate(atypeMaxCutoff)
202	+	allocate(atypeMaxCutoff(nAtypes))
203	+	endif
204
205		do i = 1, nAtypes
206		call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
#	Line 197 | Line 210 | contains
210		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
211		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
212		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
213	+	call getElementProperty(atypes, i, "is_SC", i_is_SC)
214	+	call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
215
216		do j = i, nAtypes
217
#	Line 210 | Line 225 | contains
225		call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
226		call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
227		call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
228	+	call getElementProperty(atypes, j, "is_SC", j_is_SC)
229	+	call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
230
231		if (i_is_LJ .and. j_is_LJ) then
232		iHash = ior(iHash, LJ_PAIR)
#	Line 231 | Line 248 | contains
248		iHash = ior(iHash, EAM_PAIR)
249		endif
250
251	+	if (i_is_SC .and. j_is_SC) then
252	+	iHash = ior(iHash, SC_PAIR)
253	+	endif
254	+
255		if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
256		if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
257		if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
#	Line 240 | Line 261 | contains
261		if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
262
263
264	<	InteractionMap(i,j)%InteractionHash = iHash
265	<	InteractionMap(j,i)%InteractionHash = iHash
264	>	InteractionHash(i,j) = iHash
265	>	InteractionHash(j,i) = iHash
266
267		end do
268
269		end do
249	–	end subroutine createInteractionMap
270
271	<	! 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.
272	<	subroutine createRcuts(defaultRList)
253	<	real(kind=dp), intent(in), optional :: defaultRList
254	<	integer :: iMap
255	<	integer :: map_i,map_j
256	<	real(kind=dp) :: thisRCut = 0.0_dp
257	<	real(kind=dp) :: actualCutoff = 0.0_dp
258	<	integer :: nAtypes
271	>	haveInteractionHash = .true.
272	>	end subroutine createInteractionHash
273
274	<	if(.not. allocated(InteractionMap)) return
274	>	subroutine createGtypeCutoffMap()
275
276	+	logical :: i_is_LJ
277	+	logical :: i_is_Elect
278	+	logical :: i_is_Sticky
279	+	logical :: i_is_StickyP
280	+	logical :: i_is_GB
281	+	logical :: i_is_EAM
282	+	logical :: i_is_Shape
283	+	logical :: i_is_SC
284	+	logical :: GtypeFound
285	+
286	+	integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
287	+	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
288	+	integer :: nGroupsInRow
289	+	integer :: nGroupsInCol
290	+	integer :: nGroupTypesRow,nGroupTypesCol
291	+	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
292	+	real(kind=dp) :: biggestAtypeCutoff
293	+
294	+	if (.not. haveInteractionHash) then
295	+	call createInteractionHash()
296	+	endif
297	+	#ifdef IS_MPI
298	+	nGroupsInRow = getNgroupsInRow(plan_group_row)
299	+	nGroupsInCol = getNgroupsInCol(plan_group_col)
300	+	#endif
301		nAtypes = getSize(atypes)
302	<	! If we pass a default rcut, set all atypes to that cutoff distance
303	<	if(present(defaultRList)) then
304	<	InteractionMap(:,:)%rList = defaultRList
305	<	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
306	<	haveRlist = .true.
307	<	return
308	<	end if
309	<
310	<	do map_i = 1,nAtypes
311	<	do map_j = map_i,nAtypes
312	<	iMap = InteractionMap(map_i, map_j)%InteractionHash
313	<
314	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
315	<	!            thisRCut = getLJCutOff(map_i,map_j)
316	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
302	>	! Set all of the initial cutoffs to zero.
303	>	atypeMaxCutoff = 0.0_dp
304	>	do i = 1, nAtypes
305	>	if (SimHasAtype(i)) then
306	>	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
307	>	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
308	>	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
309	>	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
310	>	call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
311	>	call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
312	>	call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
313	>	call getElementProperty(atypes, i, "is_SC", i_is_SC)
314	>
315	>	if (haveDefaultCutoffs) then
316	>	atypeMaxCutoff(i) = defaultRcut
317	>	else
318	>	if (i_is_LJ) then
319	>	thisRcut = getSigma(i) * 2.5_dp
320	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321	>	endif
322	>	if (i_is_Elect) then
323	>	thisRcut = defaultRcut
324	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325	>	endif
326	>	if (i_is_Sticky) then
327	>	thisRcut = getStickyCut(i)
328	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329	>	endif
330	>	if (i_is_StickyP) then
331	>	thisRcut = getStickyPowerCut(i)
332	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333	>	endif
334	>	if (i_is_GB) then
335	>	thisRcut = getGayBerneCut(i)
336	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337	>	endif
338	>	if (i_is_EAM) then
339	>	thisRcut = getEAMCut(i)
340	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
341	>	endif
342	>	if (i_is_Shape) then
343	>	thisRcut = getShapeCut(i)
344	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
345	>	endif
346	>	if (i_is_SC) then
347	>	thisRcut = getSCCut(i)
348	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
349	>	endif
350		endif
351	<
352	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
353	<	!            thisRCut = getElectrostaticCutOff(map_i,map_j)
282	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
351	>
352	>	if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
353	>	biggestAtypeCutoff = atypeMaxCutoff(i)
354		endif
284	–
285	–	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
286	–	!             thisRCut = getStickyCutOff(map_i,map_j)
287	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
288	–	endif
289	–
290	–	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
291	–	!              thisRCut = getStickyPowerCutOff(map_i,map_j)
292	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
293	–	endif
294	–
295	–	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
296	–	!              thisRCut = getGayberneCutOff(map_i,map_j)
297	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
298	–	endif
299	–
300	–	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
301	–	!              thisRCut = getGaybrneLJCutOff(map_i,map_j)
302	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
303	–	endif
304	–
305	–	if ( iand(iMap, EAM_PAIR).ne.0 ) then
306	–	!              thisRCut = getEAMCutOff(map_i,map_j)
307	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
308	–	endif
309	–
310	–	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
311	–	!              thisRCut = getShapeCutOff(map_i,map_j)
312	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
313	–	endif
314	–
315	–	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
316	–	!              thisRCut = getShapeLJCutOff(map_i,map_j)
317	–	if (thisRcut > actualCutoff) actualCutoff = thisRcut
318	–	endif
319	–	InteractionMap(map_i, map_j)%rList = actualCutoff
320	–	InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
321	–	end do
322	–	end do
323	–	haveRlist = .true.
324	–	end subroutine createRcuts
355
356	+	endif
357	+	enddo
358	+
359	+	istart = 1
360	+	jstart = 1
361	+	#ifdef IS_MPI
362	+	iend = nGroupsInRow
363	+	jend = nGroupsInCol
364	+	#else
365	+	iend = nGroups
366	+	jend = nGroups
367	+	#endif
368	+
369	+	!! allocate the groupToGtype and gtypeMaxCutoff here.
370	+	if(.not.allocated(groupToGtypeRow)) then
371	+	!  allocate(groupToGtype(iend))
372	+	allocate(groupToGtypeRow(iend))
373	+	else
374	+	deallocate(groupToGtypeRow)
375	+	allocate(groupToGtypeRow(iend))
376	+	endif
377	+	if(.not.allocated(groupMaxCutoffRow)) then
378	+	allocate(groupMaxCutoffRow(iend))
379	+	else
380	+	deallocate(groupMaxCutoffRow)
381	+	allocate(groupMaxCutoffRow(iend))
382	+	end if
383
384	<	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
385	<	!!$  subroutine setRlistDF( this_rlist )
386	<	!!$
387	<	!!$   real(kind=dp) :: this_rlist
388	<	!!$
389	<	!!$    rlist = this_rlist
333	<	!!$    rlistsq = rlist * rlist
334	<	!!$
335	<	!!$    haveRlist = .true.
336	<	!!$
337	<	!!$  end subroutine setRlistDF
384	>	if(.not.allocated(gtypeMaxCutoffRow)) then
385	>	allocate(gtypeMaxCutoffRow(iend))
386	>	else
387	>	deallocate(gtypeMaxCutoffRow)
388	>	allocate(gtypeMaxCutoffRow(iend))
389	>	endif
390
391
392	<	subroutine setSimVariables()
393	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
394	<	SIM_uses_LennardJones = SimUsesLennardJones()
395	<	SIM_uses_Electrostatics = SimUsesElectrostatics()
396	<	SIM_uses_Charges = SimUsesCharges()
397	<	SIM_uses_Dipoles = SimUsesDipoles()
398	<	SIM_uses_Sticky = SimUsesSticky()
399	<	SIM_uses_StickyPower = SimUsesStickyPower()
348	<	SIM_uses_GayBerne = SimUsesGayBerne()
349	<	SIM_uses_EAM = SimUsesEAM()
350	<	SIM_uses_Shapes = SimUsesShapes()
351	<	SIM_uses_FLARB = SimUsesFLARB()
352	<	SIM_uses_RF = SimUsesRF()
353	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
354	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
355	<	SIM_uses_PBC = SimUsesPBC()
392	>	#ifdef IS_MPI
393	>	! We only allocate new storage if we are in MPI because Ncol /= Nrow
394	>	if(.not.associated(groupToGtypeCol)) then
395	>	allocate(groupToGtypeCol(jend))
396	>	else
397	>	deallocate(groupToGtypeCol)
398	>	allocate(groupToGtypeCol(jend))
399	>	end if
400
401	<	haveSIMvariables = .true.
401	>	if(.not.associated(groupMaxCutoffCol)) then
402	>	allocate(groupMaxCutoffCol(jend))
403	>	else
404	>	deallocate(groupMaxCutoffCol)
405	>	allocate(groupMaxCutoffCol(jend))
406	>	end if
407	>	if(.not.associated(gtypeMaxCutoffCol)) then
408	>	allocate(gtypeMaxCutoffCol(jend))
409	>	else
410	>	deallocate(gtypeMaxCutoffCol)
411	>	allocate(gtypeMaxCutoffCol(jend))
412	>	end if
413
414	<	return
415	<	end subroutine setSimVariables
414	>	groupMaxCutoffCol = 0.0_dp
415	>	gtypeMaxCutoffCol = 0.0_dp
416
417	+	#endif
418	+	groupMaxCutoffRow = 0.0_dp
419	+	gtypeMaxCutoffRow = 0.0_dp
420	+
421	+
422	+	!! first we do a single loop over the cutoff groups to find the
423	+	!! largest cutoff for any atypes present in this group.  We also
424	+	!! create gtypes at this point.
425	+
426	+	tol = 1.0e-6_dp
427	+	nGroupTypesRow = 0
428	+	nGroupTypesCol = 0
429	+	do i = istart, iend
430	+	n_in_i = groupStartRow(i+1) - groupStartRow(i)
431	+	groupMaxCutoffRow(i) = 0.0_dp
432	+	do ia = groupStartRow(i), groupStartRow(i+1)-1
433	+	atom1 = groupListRow(ia)
434	+	#ifdef IS_MPI
435	+	me_i = atid_row(atom1)
436	+	#else
437	+	me_i = atid(atom1)
438	+	#endif
439	+	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
440	+	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
441	+	endif
442	+	enddo
443	+	if (nGroupTypesRow.eq.0) then
444	+	nGroupTypesRow = nGroupTypesRow + 1
445	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
446	+	groupToGtypeRow(i) = nGroupTypesRow
447	+	else
448	+	GtypeFound = .false.
449	+	do g = 1, nGroupTypesRow
450	+	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
451	+	groupToGtypeRow(i) = g
452	+	GtypeFound = .true.
453	+	endif
454	+	enddo
455	+	if (.not.GtypeFound) then
456	+	nGroupTypesRow = nGroupTypesRow + 1
457	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
458	+	groupToGtypeRow(i) = nGroupTypesRow
459	+	endif
460	+	endif
461	+	enddo
462	+
463	+	#ifdef IS_MPI
464	+	do j = jstart, jend
465	+	n_in_j = groupStartCol(j+1) - groupStartCol(j)
466	+	groupMaxCutoffCol(j) = 0.0_dp
467	+	do ja = groupStartCol(j), groupStartCol(j+1)-1
468	+	atom1 = groupListCol(ja)
469	+
470	+	me_j = atid_col(atom1)
471	+
472	+	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
473	+	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
474	+	endif
475	+	enddo
476	+
477	+	if (nGroupTypesCol.eq.0) then
478	+	nGroupTypesCol = nGroupTypesCol + 1
479	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
480	+	groupToGtypeCol(j) = nGroupTypesCol
481	+	else
482	+	GtypeFound = .false.
483	+	do g = 1, nGroupTypesCol
484	+	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
485	+	groupToGtypeCol(j) = g
486	+	GtypeFound = .true.
487	+	endif
488	+	enddo
489	+	if (.not.GtypeFound) then
490	+	nGroupTypesCol = nGroupTypesCol + 1
491	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
492	+	groupToGtypeCol(j) = nGroupTypesCol
493	+	endif
494	+	endif
495	+	enddo
496	+
497	+	#else
498	+	! Set pointers to information we just found
499	+	nGroupTypesCol = nGroupTypesRow
500	+	groupToGtypeCol => groupToGtypeRow
501	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
502	+	groupMaxCutoffCol => groupMaxCutoffRow
503	+	#endif
504	+
505	+	!! allocate the gtypeCutoffMap here.
506	+	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
507	+	!! then we do a double loop over all the group TYPES to find the cutoff
508	+	!! map between groups of two types
509	+	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
510	+
511	+	do i = 1, nGroupTypesRow
512	+	do j = 1, nGroupTypesCol
513	+
514	+	select case(cutoffPolicy)
515	+	case(TRADITIONAL_CUTOFF_POLICY)
516	+	thisRcut = tradRcut
517	+	case(MIX_CUTOFF_POLICY)
518	+	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
519	+	case(MAX_CUTOFF_POLICY)
520	+	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
521	+	case default
522	+	call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
523	+	return
524	+	end select
525	+	gtypeCutoffMap(i,j)%rcut = thisRcut
526	+
527	+	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
528	+
529	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
530	+
531	+	if (.not.haveSkinThickness) then
532	+	skinThickness = 1.0_dp
533	+	endif
534	+
535	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
536	+
537	+	! sanity check
538	+
539	+	if (haveDefaultCutoffs) then
540	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
541	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
542	+	endif
543	+	endif
544	+	enddo
545	+	enddo
546	+
547	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
548	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
549	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
550	+	#ifdef IS_MPI
551	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
552	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
553	+	#endif
554	+	groupMaxCutoffCol => null()
555	+	gtypeMaxCutoffCol => null()
556	+
557	+	haveGtypeCutoffMap = .true.
558	+	end subroutine createGtypeCutoffMap
559	+
560	+	subroutine setCutoffs(defRcut, defRsw)
561	+
562	+	real(kind=dp),intent(in) :: defRcut, defRsw
563	+	character(len = statusMsgSize) :: errMsg
564	+	integer :: localError
565	+
566	+	defaultRcut = defRcut
567	+	defaultRsw = defRsw
568	+
569	+	defaultDoShift = .false.
570	+	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
571	+
572	+	write(errMsg, *) &
573	+	'cutoffRadius and switchingRadius are set to the same', newline &
574	+	// tab, 'value.  OOPSE will use shifted ', newline &
575	+	// tab, 'potentials instead of switching functions.'
576	+
577	+	call handleInfo("setCutoffs", errMsg)
578	+
579	+	defaultDoShift = .true.
580	+
581	+	endif
582	+
583	+	localError = 0
584	+	call setLJDefaultCutoff( defaultRcut, defaultDoShift )
585	+	call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
586	+	call setCutoffEAM( defaultRcut )
587	+	call setCutoffSC( defaultRcut )
588	+	call set_switch(defaultRsw, defaultRcut)
589	+	call setHmatDangerousRcutValue(defaultRcut)
590	+
591	+	haveDefaultCutoffs = .true.
592	+	haveGtypeCutoffMap = .false.
593	+
594	+	end subroutine setCutoffs
595	+
596	+	subroutine cWasLame()
597	+
598	+	VisitCutoffsAfterComputing = .true.
599	+	return
600	+
601	+	end subroutine cWasLame
602	+
603	+	subroutine setCutoffPolicy(cutPolicy)
604	+
605	+	integer, intent(in) :: cutPolicy
606	+
607	+	cutoffPolicy = cutPolicy
608	+	haveCutoffPolicy = .true.
609	+	haveGtypeCutoffMap = .false.
610	+
611	+	end subroutine setCutoffPolicy
612	+
613	+	subroutine setElectrostaticMethod( thisESM )
614	+
615	+	integer, intent(in) :: thisESM
616	+
617	+	electrostaticSummationMethod = thisESM
618	+	haveElectrostaticSummationMethod = .true.
619	+
620	+	end subroutine setElectrostaticMethod
621	+
622	+	subroutine setSkinThickness( thisSkin )
623	+
624	+	real(kind=dp), intent(in) :: thisSkin
625	+
626	+	skinThickness = thisSkin
627	+	haveSkinThickness = .true.
628	+	haveGtypeCutoffMap = .false.
629	+
630	+	end subroutine setSkinThickness
631	+
632	+	subroutine setSimVariables()
633	+	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
634	+	SIM_uses_EAM = SimUsesEAM()
635	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
636	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
637	+	SIM_uses_PBC = SimUsesPBC()
638	+	SIM_uses_SC = SimUsesSC()
639	+
640	+	haveSIMvariables = .true.
641	+
642	+	return
643	+	end subroutine setSimVariables
644	+
645		subroutine doReadyCheck(error)
646		integer, intent(out) :: error
364	–
647		integer :: myStatus
648
649		error = 0
650
651	<	if (.not. haveInteractionMap) then
651	>	if (.not. haveInteractionHash) then
652	>	call createInteractionHash()
653	>	endif
654
655	<	myStatus = 0
655	>	if (.not. haveGtypeCutoffMap) then
656	>	call createGtypeCutoffMap()
657	>	endif
658
659	<	call createInteractionMap(myStatus)
660	<
661	<	if (myStatus .ne. 0) then
662	<	write(default_error, *) 'createInteractionMap failed in doForces!'
663	<	error = -1
664	<	return
379	<	endif
659	>	if (VisitCutoffsAfterComputing) then
660	>	call set_switch(largestRcut, largestRcut)
661	>	call setHmatDangerousRcutValue(largestRcut)
662	>	call setCutoffEAM(largestRcut)
663	>	call setCutoffSC(largestRcut)
664	>	VisitCutoffsAfterComputing = .false.
665		endif
666
667		if (.not. haveSIMvariables) then
668		call setSimVariables()
669		endif
670
386	–	if (.not. haveRlist) then
387	–	write(default_error, *) 'rList has not been set in doForces!'
388	–	error = -1
389	–	return
390	–	endif
391	–
671		if (.not. haveNeighborList) then
672		write(default_error, *) 'neighbor list has not been initialized in doForces!'
673		error = -1
674		return
675		end if
676	<
676	>
677		if (.not. haveSaneForceField) then
678		write(default_error, *) 'Force Field is not sane in doForces!'
679		error = -1
680		return
681		end if
682	<
682	>
683		#ifdef IS_MPI
684		if (.not. isMPISimSet()) then
685		write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
#	Line 412 | Line 691 | contains
691		end subroutine doReadyCheck
692
693
694	<	subroutine init_FF(use_RF_c, thisStat)
694	>	subroutine init_FF(thisStat)
695
417	–	logical, intent(in) :: use_RF_c
418	–
696		integer, intent(out) :: thisStat
697		integer :: my_status, nMatches
698		integer, pointer :: MatchList(:) => null()
422	–	real(kind=dp) :: rcut, rrf, rt, dielect
699
700		!! assume things are copacetic, unless they aren't
701		thisStat = 0
702
427	–	!! Fortran's version of a cast:
428	–	FF_uses_RF = use_RF_c
429	–
703		!! init_FF is called *after* all of the atom types have been
704		!! defined in atype_module using the new_atype subroutine.
705		!!
#	Line 434 | Line 707 | contains
707		!! interactions are used by the force field.
708
709		FF_uses_DirectionalAtoms = .false.
437	–	FF_uses_LennardJones = .false.
438	–	FF_uses_Electrostatics = .false.
439	–	FF_uses_Charges = .false.
710		FF_uses_Dipoles = .false.
441	–	FF_uses_Sticky = .false.
442	–	FF_uses_StickyPower = .false.
711		FF_uses_GayBerne = .false.
712		FF_uses_EAM = .false.
713	<	FF_uses_Shapes = .false.
446	<	FF_uses_FLARB = .false.
713	>	FF_uses_SC = .false.
714
715		call getMatchingElementList(atypes, "is_Directional", .true., &
716		nMatches, MatchList)
717		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
718
452	–	call getMatchingElementList(atypes, "is_LennardJones", .true., &
453	–	nMatches, MatchList)
454	–	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
455	–
456	–	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
457	–	nMatches, MatchList)
458	–	if (nMatches .gt. 0) then
459	–	FF_uses_Electrostatics = .true.
460	–	endif
461	–
462	–	call getMatchingElementList(atypes, "is_Charge", .true., &
463	–	nMatches, MatchList)
464	–	if (nMatches .gt. 0) then
465	–	FF_uses_Charges = .true.
466	–	FF_uses_Electrostatics = .true.
467	–	endif
468	–
719		call getMatchingElementList(atypes, "is_Dipole", .true., &
720		nMatches, MatchList)
721	<	if (nMatches .gt. 0) then
472	<	FF_uses_Dipoles = .true.
473	<	FF_uses_Electrostatics = .true.
474	<	FF_uses_DirectionalAtoms = .true.
475	<	endif
476	<
477	<	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
478	<	nMatches, MatchList)
479	<	if (nMatches .gt. 0) then
480	<	FF_uses_Quadrupoles = .true.
481	<	FF_uses_Electrostatics = .true.
482	<	FF_uses_DirectionalAtoms = .true.
483	<	endif
484	<
485	<	call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
486	<	MatchList)
487	<	if (nMatches .gt. 0) then
488	<	FF_uses_Sticky = .true.
489	<	FF_uses_DirectionalAtoms = .true.
490	<	endif
491	<
492	<	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
493	<	MatchList)
494	<	if (nMatches .gt. 0) then
495	<	FF_uses_StickyPower = .true.
496	<	FF_uses_DirectionalAtoms = .true.
497	<	endif
721	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
722
723		call getMatchingElementList(atypes, "is_GayBerne", .true., &
724		nMatches, MatchList)
725	<	if (nMatches .gt. 0) then
502	<	FF_uses_GayBerne = .true.
503	<	FF_uses_DirectionalAtoms = .true.
504	<	endif
725	>	if (nMatches .gt. 0) FF_uses_GayBerne = .true.
726
727		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
728		if (nMatches .gt. 0) FF_uses_EAM = .true.
729
730	<	call getMatchingElementList(atypes, "is_Shape", .true., &
731	<	nMatches, MatchList)
511	<	if (nMatches .gt. 0) then
512	<	FF_uses_Shapes = .true.
513	<	FF_uses_DirectionalAtoms = .true.
514	<	endif
730	>	call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
731	>	if (nMatches .gt. 0) FF_uses_SC = .true.
732
516	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
517	–	nMatches, MatchList)
518	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
733
520	–	!! Assume sanity (for the sake of argument)
734		haveSaneForceField = .true.
735
523	–	!! check to make sure the FF_uses_RF setting makes sense
524	–
525	–	if (FF_uses_dipoles) then
526	–	if (FF_uses_RF) then
527	–	dielect = getDielect()
528	–	call initialize_rf(dielect)
529	–	endif
530	–	else
531	–	if (FF_uses_RF) then
532	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
533	–	thisStat = -1
534	–	haveSaneForceField = .false.
535	–	return
536	–	endif
537	–	endif
538	–
539	–	!sticky module does not contain check_sticky_FF anymore
540	–	!if (FF_uses_sticky) then
541	–	!   call check_sticky_FF(my_status)
542	–	!   if (my_status /= 0) then
543	–	!      thisStat = -1
544	–	!      haveSaneForceField = .false.
545	–	!      return
546	–	!   end if
547	–	!endif
548	–
736		if (FF_uses_EAM) then
737		call init_EAM_FF(my_status)
738		if (my_status /= 0) then
#	Line 556 | Line 743 | contains
743		end if
744		endif
745
559	–	if (FF_uses_GayBerne) then
560	–	call check_gb_pair_FF(my_status)
561	–	if (my_status .ne. 0) then
562	–	thisStat = -1
563	–	haveSaneForceField = .false.
564	–	return
565	–	endif
566	–	endif
567	–
568	–	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
569	–	endif
570	–
746		if (.not. haveNeighborList) then
747		!! Create neighbor lists
748		call expandNeighborList(nLocal, my_status)
#	Line 601 | Line 776 | contains
776
777		!! Stress Tensor
778		real( kind = dp), dimension(9) :: tau
779	<	real ( kind = dp ) :: pot
779	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
780		logical ( kind = 2) :: do_pot_c, do_stress_c
781		logical :: do_pot
782		logical :: do_stress
783		logical :: in_switching_region
784		#ifdef IS_MPI
785	<	real( kind = DP ) :: pot_local
785	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
786		integer :: nAtomsInRow
787		integer :: nAtomsInCol
788		integer :: nprocs
#	Line 622 | Line 797 | contains
797		integer :: nlist
798		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
799		real( kind = DP ) :: sw, dswdr, swderiv, mf
800	+	real( kind = DP ) :: rVal
801		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
802		real(kind=dp) :: rfpot, mu_i, virial
803	+	real(kind=dp):: rCut
804		integer :: me_i, me_j, n_in_i, n_in_j
805		logical :: is_dp_i
806		integer :: neighborListSize
#	Line 631 | Line 808 | contains
808		integer :: localError
809		integer :: propPack_i, propPack_j
810		integer :: loopStart, loopEnd, loop
811	<	integer :: iMap
812	<	real(kind=dp) :: listSkin = 1.0
811	>	integer :: iHash
812	>	integer :: i1
813	>
814
815		!! initialize local variables
816
#	Line 696 | Line 874 | contains
874		! (but only on the first time through):
875		if (loop .eq. loopStart) then
876		#ifdef IS_MPI
877	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
877	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
878		update_nlist)
879		#else
880	<	call checkNeighborList(nGroups, q_group, listSkin, &
880	>	call checkNeighborList(nGroups, q_group, skinThickness, &
881		update_nlist)
882		#endif
883		endif
#	Line 750 | Line 928 | contains
928		endif
929
930		#ifdef IS_MPI
931	+	me_j = atid_col(j)
932		call get_interatomic_vector(q_group_Row(:,i), &
933		q_group_Col(:,j), d_grp, rgrpsq)
934		#else
935	+	me_j = atid(j)
936		call get_interatomic_vector(q_group(:,i), &
937		q_group(:,j), d_grp, rgrpsq)
938	<	#endif
938	>	#endif
939
940	<	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
940	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
941		if (update_nlist) then
942		nlist = nlist + 1
943
#	Line 777 | Line 957 | contains
957
958		list(nlist) = j
959		endif
960	+
961	+	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
962
963	<	if (loop .eq. PAIR_LOOP) then
964	<	vij = 0.0d0
965	<	fij(1:3) = 0.0d0
966	<	endif
967	<
968	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
969	<	in_switching_region)
970	<
971	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
972	<
973	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
974	<
975	<	atom1 = groupListRow(ia)
976	<
977	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
978	<
979	<	atom2 = groupListCol(jb)
980	<
981	<	if (skipThisPair(atom1, atom2)) cycle inner
982	<
983	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
984	<	d_atm(1:3) = d_grp(1:3)
985	<	ratmsq = rgrpsq
986	<	else
963	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
964	>	if (loop .eq. PAIR_LOOP) then
965	>	vij = 0.0_dp
966	>	fij(1) = 0.0_dp
967	>	fij(2) = 0.0_dp
968	>	fij(3) = 0.0_dp
969	>	endif
970	>
971	>	call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
972	>
973	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
974	>
975	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
976	>
977	>	atom1 = groupListRow(ia)
978	>
979	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
980	>
981	>	atom2 = groupListCol(jb)
982	>
983	>	if (skipThisPair(atom1, atom2))  cycle inner
984	>
985	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
986	>	d_atm(1) = d_grp(1)
987	>	d_atm(2) = d_grp(2)
988	>	d_atm(3) = d_grp(3)
989	>	ratmsq = rgrpsq
990	>	else
991		#ifdef IS_MPI
992	<	call get_interatomic_vector(q_Row(:,atom1), &
993	<	q_Col(:,atom2), d_atm, ratmsq)
992	>	call get_interatomic_vector(q_Row(:,atom1), &
993	>	q_Col(:,atom2), d_atm, ratmsq)
994		#else
995	<	call get_interatomic_vector(q(:,atom1), &
996	<	q(:,atom2), d_atm, ratmsq)
995	>	call get_interatomic_vector(q(:,atom1), &
996	>	q(:,atom2), d_atm, ratmsq)
997		#endif
998	<	endif
999	<
1000	<	if (loop .eq. PREPAIR_LOOP) then
998	>	endif
999	>
1000	>	if (loop .eq. PREPAIR_LOOP) then
1001		#ifdef IS_MPI
1002	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1003	<	rgrpsq, d_grp, do_pot, do_stress, &
1004	<	eFrame, A, f, t, pot_local)
1002	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1003	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1004	>	eFrame, A, f, t, pot_local)
1005		#else
1006	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1007	<	rgrpsq, d_grp, do_pot, do_stress, &
1008	<	eFrame, A, f, t, pot)
1009	<	#endif
1010	<	else
1006	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1007	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1008	>	eFrame, A, f, t, pot)
1009	>	#endif
1010	>	else
1011		#ifdef IS_MPI
1012	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1013	<	do_pot, &
1014	<	eFrame, A, f, t, pot_local, vpair, fpair)
1012	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1013	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1014	>	fpair, d_grp, rgrp, rCut)
1015		#else
1016	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1017	<	do_pot,  &
1018	<	eFrame, A, f, t, pot, vpair, fpair)
1016	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1017	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1018	>	d_grp, rgrp, rCut)
1019		#endif
1020	+	vij = vij + vpair
1021	+	fij(1) = fij(1) + fpair(1)
1022	+	fij(2) = fij(2) + fpair(2)
1023	+	fij(3) = fij(3) + fpair(3)
1024	+	endif
1025	+	enddo inner
1026	+	enddo
1027
1028	<	vij = vij + vpair
1029	<	fij(1:3) = fij(1:3) + fpair(1:3)
1030	<	endif
1031	<	enddo inner
1032	<	enddo
1033	<
1034	<	if (loop .eq. PAIR_LOOP) then
1035	<	if (in_switching_region) then
1036	<	swderiv = vij*dswdr/rgrp
1037	<	fij(1) = fij(1) + swderiv*d_grp(1)
845	<	fij(2) = fij(2) + swderiv*d_grp(2)
846	<	fij(3) = fij(3) + swderiv*d_grp(3)
847	<
848	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
849	<	atom1=groupListRow(ia)
850	<	mf = mfactRow(atom1)
1028	>	if (loop .eq. PAIR_LOOP) then
1029	>	if (in_switching_region) then
1030	>	swderiv = vij*dswdr/rgrp
1031	>	fij(1) = fij(1) + swderiv*d_grp(1)
1032	>	fij(2) = fij(2) + swderiv*d_grp(2)
1033	>	fij(3) = fij(3) + swderiv*d_grp(3)
1034	>
1035	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1036	>	atom1=groupListRow(ia)
1037	>	mf = mfactRow(atom1)
1038		#ifdef IS_MPI
1039	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1040	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1041	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1039	>	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1040	>	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1041	>	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1042		#else
1043	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1044	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1045	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1043	>	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1044	>	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1045	>	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1046		#endif
1047	<	enddo
1048	<
1049	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1050	<	atom2=groupListCol(jb)
1051	<	mf = mfactCol(atom2)
1047	>	enddo
1048	>
1049	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1050	>	atom2=groupListCol(jb)
1051	>	mf = mfactCol(atom2)
1052		#ifdef IS_MPI
1053	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1054	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1055	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1053	>	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1054	>	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1055	>	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1056		#else
1057	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1058	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1059	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1057	>	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1058	>	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1059	>	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1060		#endif
1061	<	enddo
1062	<	endif
1061	>	enddo
1062	>	endif
1063
1064	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1064	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1065	>	endif
1066		endif
1067	<	end if
1067	>	endif
1068		enddo
1069	+
1070		enddo outer
1071
1072		if (update_nlist) then
#	Line 937 | Line 1126 | contains
1126
1127		if (do_pot) then
1128		! scatter/gather pot_row into the members of my column
1129	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1130	<
1129	>	do i = 1,LR_POT_TYPES
1130	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1131	>	end do
1132		! scatter/gather pot_local into all other procs
1133		! add resultant to get total pot
1134		do i = 1, nlocal
1135	<	pot_local = pot_local + pot_Temp(i)
1135	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1136	>	+ pot_Temp(1:LR_POT_TYPES,i)
1137		enddo
1138
1139		pot_Temp = 0.0_DP
1140	<
1141	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1140	>	do i = 1,LR_POT_TYPES
1141	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1142	>	end do
1143		do i = 1, nlocal
1144	<	pot_local = pot_local + pot_Temp(i)
1144	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1145	>	+ pot_Temp(1:LR_POT_TYPES,i)
1146		enddo
1147
1148		endif
1149		#endif
1150
1151	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1151	>	if (SIM_requires_postpair_calc) then
1152	>	do i = 1, nlocal
1153	>
1154	>	! we loop only over the local atoms, so we don't need row and column
1155	>	! lookups for the types
1156	>
1157	>	me_i = atid(i)
1158	>
1159	>	! is the atom electrostatic?  See if it would have an
1160	>	! electrostatic interaction with itself
1161	>	iHash = InteractionHash(me_i,me_i)
1162
1163	<	if (FF_uses_RF .and. SIM_uses_RF) then
961	<
1163	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1164		#ifdef IS_MPI
1165	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1166	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
965	<	do i = 1,nlocal
966	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
967	<	end do
968	<	#endif
969	<
970	<	do i = 1, nLocal
971	<
972	<	rfpot = 0.0_DP
973	<	#ifdef IS_MPI
974	<	me_i = atid_row(i)
1165	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1166	>	t, do_pot)
1167		#else
1168	<	me_i = atid(i)
1168	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1169	>	t, do_pot)
1170		#endif
1171	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1171	>	endif
1172	>
1173	>
1174	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1175
1176	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1177	<
1178	<	mu_i = getDipoleMoment(me_i)
1179	<
1180	<	!! The reaction field needs to include a self contribution
1181	<	!! to the field:
1182	<	call accumulate_self_rf(i, mu_i, eFrame)
1183	<	!! Get the reaction field contribution to the
1184	<	!! potential and torques:
1185	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1176	>	! loop over the excludes to accumulate RF stuff we've
1177	>	! left out of the normal pair loop
1178	>
1179	>	do i1 = 1, nSkipsForAtom(i)
1180	>	j = skipsForAtom(i, i1)
1181	>
1182	>	! prevent overcounting of the skips
1183	>	if (i.lt.j) then
1184	>	call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1185	>	rVal = sqrt(ratmsq)
1186	>	call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1187		#ifdef IS_MPI
1188	<	pot_local = pot_local + rfpot
1188	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1189	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1190		#else
1191	<	pot = pot + rfpot
1192	<
1191	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1192	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1193		#endif
1194	<	endif
1195	<	enddo
1196	<	endif
1194	>	endif
1195	>	enddo
1196	>	endif
1197	>	enddo
1198		endif
1199	<
1001	<
1199	>
1200		#ifdef IS_MPI
1201	<
1201	>
1202		if (do_pot) then
1203	<	pot = pot + pot_local
1204	<	!! we assume the c code will do the allreduce to get the total potential
1205	<	!! we could do it right here if we needed to...
1203	>	#ifdef SINGLE_PRECISION
1204	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1205	>	mpi_comm_world,mpi_err)
1206	>	#else
1207	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1208	>	mpi_comm_world,mpi_err)
1209	>	#endif
1210		endif
1211	<
1211	>
1212		if (do_stress) then
1213	+	#ifdef SINGLE_PRECISION
1214	+	call mpi_allreduce(tau_Temp, tau, 9,mpi_real,mpi_sum, &
1215	+	mpi_comm_world,mpi_err)
1216	+	call mpi_allreduce(virial_Temp, virial,1,mpi_real,mpi_sum, &
1217	+	mpi_comm_world,mpi_err)
1218	+	#else
1219		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1220		mpi_comm_world,mpi_err)
1221		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1222		mpi_comm_world,mpi_err)
1223	+	#endif
1224		endif
1225	<
1225	>
1226		#else
1227	<
1227	>
1228		if (do_stress) then
1229		tau = tau_Temp
1230		virial = virial_Temp
1231		endif
1232	<
1232	>
1233		#endif
1234	<
1234	>
1235		end subroutine do_force_loop
1236
1237		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1238	<	eFrame, A, f, t, pot, vpair, fpair)
1238	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1239
1240	<	real( kind = dp ) :: pot, vpair, sw
1240	>	real( kind = dp ) :: vpair, sw
1241	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1242		real( kind = dp ), dimension(3) :: fpair
1243		real( kind = dp ), dimension(nLocal)   :: mfact
1244		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1039 | Line 1249 | contains
1249		logical, intent(inout) :: do_pot
1250		integer, intent(in) :: i, j
1251		real ( kind = dp ), intent(inout) :: rijsq
1252	<	real ( kind = dp )                :: r
1252	>	real ( kind = dp ), intent(inout) :: r_grp
1253		real ( kind = dp ), intent(inout) :: d(3)
1254	<	real ( kind = dp ) :: ebalance
1254	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1255	>	real ( kind = dp ), intent(inout) :: rCut
1256	>	real ( kind = dp ) :: r
1257	>	real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1258		integer :: me_i, me_j
1259	+	integer :: k
1260
1261	<	integer :: iMap
1261	>	integer :: iHash
1262
1263		r = sqrt(rijsq)
1264	<	vpair = 0.0d0
1265	<	fpair(1:3) = 0.0d0
1264	>
1265	>	vpair = 0.0_dp
1266	>	fpair(1:3) = 0.0_dp
1267
1268		#ifdef IS_MPI
1269		me_i = atid_row(i)
#	Line 1058 | Line 1273 | contains
1273		me_j = atid(j)
1274		#endif
1275
1276	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1277	<
1278	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1279	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1276	>	iHash = InteractionHash(me_i, me_j)
1277	>
1278	>	if ( iand(iHash, LJ_PAIR).ne.0 ) then
1279	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1280	>	pot(VDW_POT), f, do_pot)
1281		endif
1282	<
1283	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1284	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1285	<	pot, eFrame, f, t, do_pot)
1070	<
1071	<	if (FF_uses_RF .and. SIM_uses_RF) then
1072	<
1073	<	! CHECK ME (RF needs to know about all electrostatic types)
1074	<	call accumulate_rf(i, j, r, eFrame, sw)
1075	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1076	<	endif
1077	<
1282	>
1283	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1284	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1285	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1286		endif
1287	<
1288	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1287	>
1288	>	if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1289		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1290	<	pot, A, f, t, do_pot)
1290	>	pot(HB_POT), A, f, t, do_pot)
1291		endif
1292	<
1293	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1292	>
1293	>	if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1294		call do_sticky_power_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, GAYBERNE_PAIR).ne.0 ) then
1297	>
1298	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1299		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1300	<	pot, A, f, t, do_pot)
1300	>	pot(VDW_POT), A, f, t, do_pot)
1301		endif
1302
1303	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1304	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1305	<	!           pot, A, f, t, do_pot)
1303	>	if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1304	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1305	>	pot(VDW_POT), A, f, t, do_pot)
1306		endif
1307	<
1308	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1309	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1310	<	do_pot)
1307	>
1308	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1309	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1310	>	pot(METALLIC_POT), f, do_pot)
1311		endif
1312	<
1313	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1312	>
1313	>	if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1314		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1315	<	pot, A, f, t, do_pot)
1315	>	pot(VDW_POT), A, f, t, do_pot)
1316		endif
1317	<
1318	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1317	>
1318	>	if ( iand(iHash, SHAPE_LJ).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	<
1322	>
1323	>	if ( iand(iHash, SC_PAIR).ne.0 ) then
1324	>	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1325	>	pot(METALLIC_POT), f, do_pot)
1326	>	endif
1327	>
1328		end subroutine do_pair
1329
1330	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1330	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1331		do_pot, do_stress, eFrame, A, f, t, pot)
1332
1333	<	real( kind = dp ) :: pot, sw
1333	>	real( kind = dp ) :: sw
1334	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1335		real( kind = dp ), dimension(9,nLocal) :: eFrame
1336		real (kind=dp), dimension(9,nLocal) :: A
1337		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1125 | Line 1339 | contains
1339
1340		logical, intent(inout) :: do_pot, do_stress
1341		integer, intent(in) :: i, j
1342	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1342	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1343		real ( kind = dp )                :: r, rc
1344		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1345
1346	<	integer :: me_i, me_j, iMap
1346	>	integer :: me_i, me_j, iHash
1347
1348	+	r = sqrt(rijsq)
1349	+
1350		#ifdef IS_MPI
1351		me_i = atid_row(i)
1352		me_j = atid_col(j)
#	Line 1139 | Line 1355 | contains
1355		me_j = atid(j)
1356		#endif
1357
1358	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1358	>	iHash = InteractionHash(me_i, me_j)
1359
1360	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1361	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1360	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1361	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1362		endif
1363	+
1364	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1365	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1366	+	endif
1367
1368		end subroutine do_prepair
1369
1370
1371		subroutine do_preforce(nlocal,pot)
1372		integer :: nlocal
1373	<	real( kind = dp ) :: pot
1373	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1374
1375		if (FF_uses_EAM .and. SIM_uses_EAM) then
1376	<	call calc_EAM_preforce_Frho(nlocal,pot)
1376	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1377		endif
1378	<
1379	<
1378	>	if (FF_uses_SC .and. SIM_uses_SC) then
1379	>	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1380	>	endif
1381		end subroutine do_preforce
1382
1383
#	Line 1168 | Line 1389 | contains
1389		real( kind = dp ) :: d(3), scaled(3)
1390		integer i
1391
1392	<	d(1:3) = q_j(1:3) - q_i(1:3)
1392	>	d(1) = q_j(1) - q_i(1)
1393	>	d(2) = q_j(2) - q_i(2)
1394	>	d(3) = q_j(3) - q_i(3)
1395
1396		! Wrap back into periodic box if necessary
1397		if ( SIM_uses_PBC ) then
1398
1399		if( .not.boxIsOrthorhombic ) then
1400		! calc the scaled coordinates.
1401	<
1179	<	scaled = matmul(HmatInv, d)
1401	>	! scaled = matmul(HmatInv, d)
1402
1403	+	scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1404	+	scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1405	+	scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1406	+
1407		! wrap the scaled coordinates
1408
1409	<	scaled = scaled  - anint(scaled)
1409	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1410	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1411	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1412
1185	–
1413		! calc the wrapped real coordinates from the wrapped scaled
1414		! coordinates
1415	+	! d = matmul(Hmat,scaled)
1416	+	d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1417	+	d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1418	+	d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1419
1189	–	d = matmul(Hmat,scaled)
1190	–
1420		else
1421		! calc the scaled coordinates.
1422
1423	<	do i = 1, 3
1424	<	scaled(i) = d(i) * HmatInv(i,i)
1423	>	scaled(1) = d(1) * HmatInv(1,1)
1424	>	scaled(2) = d(2) * HmatInv(2,2)
1425	>	scaled(3) = d(3) * HmatInv(3,3)
1426	>
1427	>	! wrap the scaled coordinates
1428	>
1429	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1430	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1431	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1432
1433	<	! wrap the scaled coordinates
1433	>	! calc the wrapped real coordinates from the wrapped scaled
1434	>	! coordinates
1435
1436	<	scaled(i) = scaled(i) - anint(scaled(i))
1436	>	d(1) = scaled(1)*Hmat(1,1)
1437	>	d(2) = scaled(2)*Hmat(2,2)
1438	>	d(3) = scaled(3)*Hmat(3,3)
1439
1201	–	! calc the wrapped real coordinates from the wrapped scaled
1202	–	! coordinates
1203	–
1204	–	d(i) = scaled(i)*Hmat(i,i)
1205	–	enddo
1440		endif
1441
1442		endif
1443
1444	<	r_sq = dot_product(d,d)
1444	>	r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1445
1446		end subroutine get_interatomic_vector
1447
#	Line 1239 | Line 1473 | contains
1473		pot_Col = 0.0_dp
1474		pot_Temp = 0.0_dp
1475
1242	–	rf_Row = 0.0_dp
1243	–	rf_Col = 0.0_dp
1244	–	rf_Temp = 0.0_dp
1245	–
1476		#endif
1477
1478		if (FF_uses_EAM .and. SIM_uses_EAM) then
1479		call clean_EAM()
1480		endif
1481
1252	–	rf = 0.0_dp
1482		tau_Temp = 0.0_dp
1483		virial_Temp = 0.0_dp
1484		end subroutine zero_work_arrays
#	Line 1338 | Line 1567 | contains
1567
1568		function FF_UsesDirectionalAtoms() result(doesit)
1569		logical :: doesit
1570	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1342	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1343	<	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1570	>	doesit = FF_uses_DirectionalAtoms
1571		end function FF_UsesDirectionalAtoms
1572
1573		function FF_RequiresPrepairCalc() result(doesit)
1574		logical :: doesit
1575	<	doesit = FF_uses_EAM
1575	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1576	>	.or. FF_uses_MEAM
1577		end function FF_RequiresPrepairCalc
1578
1351	–	function FF_RequiresPostpairCalc() result(doesit)
1352	–	logical :: doesit
1353	–	doesit = FF_uses_RF
1354	–	end function FF_RequiresPostpairCalc
1355	–
1579		#ifdef PROFILE
1580		function getforcetime() result(totalforcetime)
1581		real(kind=dp) :: totalforcetime

Diff Legend

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