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