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

Diff Legend

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