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 2461 by gezelter, Mon Nov 21 22:59:02 2005 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.69 2005-11-21 22:58:35 gezelter Exp $, $Date: 2005-11-21 22:58:35 $, $Name: not supported by cvs2svn $, $Revision: 1.69 $
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	>	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 :: GtypeFound
285	>
286	>	integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
287	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
288	>	integer :: nGroupsInRow
289	>	integer :: nGroupsInCol
290	>	integer :: nGroupTypesRow,nGroupTypesCol
291	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
292	>	real(kind=dp) :: biggestAtypeCutoff
293
294	+	if (.not. haveInteractionHash) then
295	+	call createInteractionHash()
296	+	endif
297	+	#ifdef IS_MPI
298	+	nGroupsInRow = getNgroupsInRow(plan_group_row)
299	+	nGroupsInCol = getNgroupsInCol(plan_group_col)
300	+	#endif
301		nAtypes = getSize(atypes)
302	<	! If we pass a default rcut, set all atypes to that cutoff distance
303	<	if(present(defaultRList)) then
304	<	InteractionMap(:,:)%rList = defaultRList
305	<	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
306	<	haveRlist = .true.
307	<	return
308	<	end if
309	<
310	<	do map_i = 1,nAtypes
311	<	do map_j = map_i,nAtypes
312	<	iMap = InteractionMap(map_i, map_j)%InteractionHash
302	>	! Set all of the initial cutoffs to zero.
303	>	atypeMaxCutoff = 0.0_dp
304	>	do i = 1, nAtypes
305	>	if (SimHasAtype(i)) then
306	>	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
307	>	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
308	>	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
309	>	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
310	>	call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
311	>	call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
312	>	call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
313
314	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
315	<	!            thisRCut = getLJCutOff(map_i,map_j)
316	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
314	>
315	>	if (haveDefaultCutoffs) then
316	>	atypeMaxCutoff(i) = defaultRcut
317	>	else
318	>	if (i_is_LJ) then
319	>	thisRcut = getSigma(i) * 2.5_dp
320	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321	>	endif
322	>	if (i_is_Elect) then
323	>	thisRcut = defaultRcut
324	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325	>	endif
326	>	if (i_is_Sticky) then
327	>	thisRcut = getStickyCut(i)
328	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329	>	endif
330	>	if (i_is_StickyP) then
331	>	thisRcut = getStickyPowerCut(i)
332	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333	>	endif
334	>	if (i_is_GB) then
335	>	thisRcut = getGayBerneCut(i)
336	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337	>	endif
338	>	if (i_is_EAM) then
339	>	thisRcut = getEAMCut(i)
340	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
341	>	endif
342	>	if (i_is_Shape) then
343	>	thisRcut = getShapeCut(i)
344	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
345	>	endif
346		endif
347	<
348	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
349	<	!            thisRCut = getElectrostaticCutOff(map_i,map_j)
282	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
347	>
348	>	if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
349	>	biggestAtypeCutoff = atypeMaxCutoff(i)
350		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
351
352	+	endif
353	+	enddo
354	+
355	+	istart = 1
356	+	jstart = 1
357	+	#ifdef IS_MPI
358	+	iend = nGroupsInRow
359	+	jend = nGroupsInCol
360	+	#else
361	+	iend = nGroups
362	+	jend = nGroups
363	+	#endif
364	+
365	+	!! allocate the groupToGtype and gtypeMaxCutoff here.
366	+	if(.not.allocated(groupToGtypeRow)) then
367	+	!  allocate(groupToGtype(iend))
368	+	allocate(groupToGtypeRow(iend))
369	+	else
370	+	deallocate(groupToGtypeRow)
371	+	allocate(groupToGtypeRow(iend))
372	+	endif
373	+	if(.not.allocated(groupMaxCutoffRow)) then
374	+	allocate(groupMaxCutoffRow(iend))
375	+	else
376	+	deallocate(groupMaxCutoffRow)
377	+	allocate(groupMaxCutoffRow(iend))
378	+	end if
379
380	<	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
381	<	!!$  subroutine setRlistDF( this_rlist )
382	<	!!$
383	<	!!$   real(kind=dp) :: this_rlist
384	<	!!$
385	<	!!$    rlist = this_rlist
333	<	!!$    rlistsq = rlist * rlist
334	<	!!$
335	<	!!$    haveRlist = .true.
336	<	!!$
337	<	!!$  end subroutine setRlistDF
380	>	if(.not.allocated(gtypeMaxCutoffRow)) then
381	>	allocate(gtypeMaxCutoffRow(iend))
382	>	else
383	>	deallocate(gtypeMaxCutoffRow)
384	>	allocate(gtypeMaxCutoffRow(iend))
385	>	endif
386
387
388	<	subroutine setSimVariables()
389	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
390	<	SIM_uses_LennardJones = SimUsesLennardJones()
391	<	SIM_uses_Electrostatics = SimUsesElectrostatics()
392	<	SIM_uses_Charges = SimUsesCharges()
393	<	SIM_uses_Dipoles = SimUsesDipoles()
394	<	SIM_uses_Sticky = SimUsesSticky()
395	<	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()
388	>	#ifdef IS_MPI
389	>	! We only allocate new storage if we are in MPI because Ncol /= Nrow
390	>	if(.not.associated(groupToGtypeCol)) then
391	>	allocate(groupToGtypeCol(jend))
392	>	else
393	>	deallocate(groupToGtypeCol)
394	>	allocate(groupToGtypeCol(jend))
395	>	end if
396
397	<	haveSIMvariables = .true.
397	>	if(.not.associated(groupToGtypeCol)) then
398	>	allocate(groupToGtypeCol(jend))
399	>	else
400	>	deallocate(groupToGtypeCol)
401	>	allocate(groupToGtypeCol(jend))
402	>	end if
403	>	if(.not.associated(gtypeMaxCutoffCol)) then
404	>	allocate(gtypeMaxCutoffCol(jend))
405	>	else
406	>	deallocate(gtypeMaxCutoffCol)
407	>	allocate(gtypeMaxCutoffCol(jend))
408	>	end if
409
410	<	return
411	<	end subroutine setSimVariables
410	>	groupMaxCutoffCol = 0.0_dp
411	>	gtypeMaxCutoffCol = 0.0_dp
412
413	+	#endif
414	+	groupMaxCutoffRow = 0.0_dp
415	+	gtypeMaxCutoffRow = 0.0_dp
416	+
417	+
418	+	!! first we do a single loop over the cutoff groups to find the
419	+	!! largest cutoff for any atypes present in this group.  We also
420	+	!! create gtypes at this point.
421	+
422	+	tol = 1.0d-6
423	+	nGroupTypesRow = 0
424	+
425	+	do i = istart, iend
426	+	n_in_i = groupStartRow(i+1) - groupStartRow(i)
427	+	groupMaxCutoffRow(i) = 0.0_dp
428	+	do ia = groupStartRow(i), groupStartRow(i+1)-1
429	+	atom1 = groupListRow(ia)
430	+	#ifdef IS_MPI
431	+	me_i = atid_row(atom1)
432	+	#else
433	+	me_i = atid(atom1)
434	+	#endif
435	+	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
436	+	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
437	+	endif
438	+	enddo
439	+
440	+	if (nGroupTypesRow.eq.0) then
441	+	nGroupTypesRow = nGroupTypesRow + 1
442	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
443	+	groupToGtypeRow(i) = nGroupTypesRow
444	+	else
445	+	GtypeFound = .false.
446	+	do g = 1, nGroupTypesRow
447	+	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
448	+	groupToGtypeRow(i) = g
449	+	GtypeFound = .true.
450	+	endif
451	+	enddo
452	+	if (.not.GtypeFound) then
453	+	nGroupTypesRow = nGroupTypesRow + 1
454	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
455	+	groupToGtypeRow(i) = nGroupTypesRow
456	+	endif
457	+	endif
458	+	enddo
459	+
460	+	#ifdef IS_MPI
461	+	do j = jstart, jend
462	+	n_in_j = groupStartCol(j+1) - groupStartCol(j)
463	+	groupMaxCutoffCol(j) = 0.0_dp
464	+	do ja = groupStartCol(j), groupStartCol(j+1)-1
465	+	atom1 = groupListCol(ja)
466	+
467	+	me_j = atid_col(atom1)
468	+
469	+	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
470	+	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
471	+	endif
472	+	enddo
473	+
474	+	if (nGroupTypesCol.eq.0) then
475	+	nGroupTypesCol = nGroupTypesCol + 1
476	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
477	+	groupToGtypeCol(j) = nGroupTypesCol
478	+	else
479	+	GtypeFound = .false.
480	+	do g = 1, nGroupTypesCol
481	+	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
482	+	groupToGtypeCol(j) = g
483	+	GtypeFound = .true.
484	+	endif
485	+	enddo
486	+	if (.not.GtypeFound) then
487	+	nGroupTypesCol = nGroupTypesCol + 1
488	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
489	+	groupToGtypeCol(j) = nGroupTypesCol
490	+	endif
491	+	endif
492	+	enddo
493	+
494	+	#else
495	+	! Set pointers to information we just found
496	+	nGroupTypesCol = nGroupTypesRow
497	+	groupToGtypeCol => groupToGtypeRow
498	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
499	+	groupMaxCutoffCol => groupMaxCutoffRow
500	+	#endif
501	+
502	+	!! allocate the gtypeCutoffMap here.
503	+	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
504	+	!! then we do a double loop over all the group TYPES to find the cutoff
505	+	!! map between groups of two types
506	+	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507	+
508	+	do i = 1, nGroupTypesRow
509	+	do j = 1, nGroupTypesCol
510	+
511	+	select case(cutoffPolicy)
512	+	case(TRADITIONAL_CUTOFF_POLICY)
513	+	thisRcut = tradRcut
514	+	case(MIX_CUTOFF_POLICY)
515	+	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516	+	case(MAX_CUTOFF_POLICY)
517	+	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
518	+	case default
519	+	call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520	+	return
521	+	end select
522	+	gtypeCutoffMap(i,j)%rcut = thisRcut
523	+
524	+	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
525	+
526	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
527	+
528	+	if (.not.haveSkinThickness) then
529	+	skinThickness = 1.0_dp
530	+	endif
531	+
532	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
533	+
534	+	! sanity check
535	+
536	+	if (haveDefaultCutoffs) then
537	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
538	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
539	+	endif
540	+	endif
541	+	enddo
542	+	enddo
543	+
544	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
545	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
546	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
547	+	#ifdef IS_MPI
548	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
549	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
550	+	#endif
551	+	groupMaxCutoffCol => null()
552	+	gtypeMaxCutoffCol => null()
553	+
554	+	haveGtypeCutoffMap = .true.
555	+	end subroutine createGtypeCutoffMap
556	+
557	+	subroutine setCutoffs(defRcut, defRsw)
558	+
559	+	real(kind=dp),intent(in) :: defRcut, defRsw
560	+	character(len = statusMsgSize) :: errMsg
561	+	integer :: localError
562	+
563	+	defaultRcut = defRcut
564	+	defaultRsw = defRsw
565	+
566	+	defaultDoShift = .false.
567	+	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
568	+
569	+	write(errMsg, *) &
570	+	'cutoffRadius and switchingRadius are set to the same', newline &
571	+	// tab, 'value.  OOPSE will use shifted ', newline &
572	+	// tab, 'potentials instead of switching functions.'
573	+
574	+	call handleInfo("setCutoffs", errMsg)
575	+
576	+	defaultDoShift = .true.
577	+
578	+	endif
579	+
580	+	localError = 0
581	+	call setLJDefaultCutoff( defaultRcut, defaultDoShift )
582	+	call setCutoffEAM( defaultRcut, localError)
583	+	if (localError /= 0) then
584	+	write(errMsg, *) 'An error has occured in setting the EAM cutoff'
585	+	call handleError("setCutoffs", errMsg)
586	+	end if
587	+	call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
588	+
589	+	haveDefaultCutoffs = .true.
590	+	end subroutine setCutoffs
591	+
592	+	subroutine cWasLame()
593	+
594	+	VisitCutoffsAfterComputing = .true.
595	+	return
596	+
597	+	end subroutine cWasLame
598	+
599	+	subroutine setCutoffPolicy(cutPolicy)
600	+
601	+	integer, intent(in) :: cutPolicy
602	+
603	+	cutoffPolicy = cutPolicy
604	+	haveCutoffPolicy = .true.
605	+
606	+	call createGtypeCutoffMap()
607	+
608	+	end subroutine setCutoffPolicy
609	+
610	+	subroutine setElectrostaticMethod( thisESM )
611	+
612	+	integer, intent(in) :: thisESM
613	+
614	+	electrostaticSummationMethod = thisESM
615	+	haveElectrostaticSummationMethod = .true.
616	+
617	+	end subroutine setElectrostaticMethod
618	+
619	+	subroutine setSkinThickness( thisSkin )
620	+
621	+	real(kind=dp), intent(in) :: thisSkin
622	+
623	+	skinThickness = thisSkin
624	+	haveSkinThickness = .true.
625	+
626	+	call createGtypeCutoffMap()
627	+
628	+	end subroutine setSkinThickness
629	+
630	+	subroutine setSimVariables()
631	+	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
632	+	SIM_uses_EAM = SimUsesEAM()
633	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
634	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
635	+	SIM_uses_PBC = SimUsesPBC()
636	+
637	+	haveSIMvariables = .true.
638	+
639	+	return
640	+	end subroutine setSimVariables
641	+
642		subroutine doReadyCheck(error)
643		integer, intent(out) :: error
644
#	Line 366 | Line 646 | contains
646
647		error = 0
648
649	<	if (.not. haveInteractionMap) then
649	>	if (.not. haveInteractionHash) then
650	>	call createInteractionHash()
651	>	endif
652
653	<	myStatus = 0
653	>	if (.not. haveGtypeCutoffMap) then
654	>	call createGtypeCutoffMap()
655	>	endif
656
373	–	call createInteractionMap(myStatus)
657
658	<	if (myStatus .ne. 0) then
659	<	write(default_error, *) 'createInteractionMap failed in doForces!'
377	<	error = -1
378	<	return
379	<	endif
658	>	if (VisitCutoffsAfterComputing) then
659	>	call set_switch(GROUP_SWITCH, largestRcut, largestRcut)
660		endif
661
662	+
663		if (.not. haveSIMvariables) then
664		call setSimVariables()
665		endif
666
667	<	if (.not. haveRlist) then
668	<	write(default_error, *) 'rList has not been set in doForces!'
669	<	error = -1
670	<	return
671	<	endif
667	>	!  if (.not. haveRlist) then
668	>	!     write(default_error, *) 'rList has not been set in doForces!'
669	>	!     error = -1
670	>	!     return
671	>	!  endif
672
673		if (.not. haveNeighborList) then
674		write(default_error, *) 'neighbor list has not been initialized in doForces!'
#	Line 412 | Line 693 | contains
693		end subroutine doReadyCheck
694
695
696	<	subroutine init_FF(use_RF_c, thisStat)
696	>	subroutine init_FF(thisStat)
697
417	–	logical, intent(in) :: use_RF_c
418	–
698		integer, intent(out) :: thisStat
699		integer :: my_status, nMatches
700		integer, pointer :: MatchList(:) => null()
422	–	real(kind=dp) :: rcut, rrf, rt, dielect
701
702		!! assume things are copacetic, unless they aren't
703		thisStat = 0
704
427	–	!! Fortran's version of a cast:
428	–	FF_uses_RF = use_RF_c
429	–
705		!! init_FF is called *after* all of the atom types have been
706		!! defined in atype_module using the new_atype subroutine.
707		!!
#	Line 434 | Line 709 | contains
709		!! interactions are used by the force field.
710
711		FF_uses_DirectionalAtoms = .false.
437	–	FF_uses_LennardJones = .false.
438	–	FF_uses_Electrostatics = .false.
439	–	FF_uses_Charges = .false.
712		FF_uses_Dipoles = .false.
441	–	FF_uses_Sticky = .false.
442	–	FF_uses_StickyPower = .false.
713		FF_uses_GayBerne = .false.
714		FF_uses_EAM = .false.
445	–	FF_uses_Shapes = .false.
446	–	FF_uses_FLARB = .false.
715
716		call getMatchingElementList(atypes, "is_Directional", .true., &
717		nMatches, MatchList)
718		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
719
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	–
720		call getMatchingElementList(atypes, "is_Dipole", .true., &
721		nMatches, MatchList)
722	<	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
722	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
723
724		call getMatchingElementList(atypes, "is_GayBerne", .true., &
725		nMatches, MatchList)
726	<	if (nMatches .gt. 0) then
502	<	FF_uses_GayBerne = .true.
503	<	FF_uses_DirectionalAtoms = .true.
504	<	endif
726	>	if (nMatches .gt. 0) FF_uses_GayBerne = .true.
727
728		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
729		if (nMatches .gt. 0) FF_uses_EAM = .true.
730
509	–	call getMatchingElementList(atypes, "is_Shape", .true., &
510	–	nMatches, MatchList)
511	–	if (nMatches .gt. 0) then
512	–	FF_uses_Shapes = .true.
513	–	FF_uses_DirectionalAtoms = .true.
514	–	endif
731
516	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
517	–	nMatches, MatchList)
518	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
519	–
520	–	!! Assume sanity (for the sake of argument)
732		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
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
733
734		if (FF_uses_EAM) then
735		call init_EAM_FF(my_status)
#	Line 556 | Line 741 | contains
741		end if
742		endif
743
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	–
744		if (.not. haveNeighborList) then
745		!! Create neighbor lists
746		call expandNeighborList(nLocal, my_status)
#	Line 601 | Line 774 | contains
774
775		!! Stress Tensor
776		real( kind = dp), dimension(9) :: tau
777	<	real ( kind = dp ) :: pot
777	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
778		logical ( kind = 2) :: do_pot_c, do_stress_c
779		logical :: do_pot
780		logical :: do_stress
781		logical :: in_switching_region
782		#ifdef IS_MPI
783	<	real( kind = DP ) :: pot_local
783	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
784		integer :: nAtomsInRow
785		integer :: nAtomsInCol
786		integer :: nprocs
#	Line 622 | Line 795 | contains
795		integer :: nlist
796		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
797		real( kind = DP ) :: sw, dswdr, swderiv, mf
798	+	real( kind = DP ) :: rVal
799		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
800		real(kind=dp) :: rfpot, mu_i, virial
801	+	real(kind=dp):: rCut
802		integer :: me_i, me_j, n_in_i, n_in_j
803		logical :: is_dp_i
804		integer :: neighborListSize
#	Line 631 | Line 806 | contains
806		integer :: localError
807		integer :: propPack_i, propPack_j
808		integer :: loopStart, loopEnd, loop
809	<	integer :: iMap
810	<	real(kind=dp) :: listSkin = 1.0
809	>	integer :: iHash
810	>	integer :: i1
811	>
812
813		!! initialize local variables
814
#	Line 696 | Line 872 | contains
872		! (but only on the first time through):
873		if (loop .eq. loopStart) then
874		#ifdef IS_MPI
875	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
875	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
876		update_nlist)
877		#else
878	<	call checkNeighborList(nGroups, q_group, listSkin, &
878	>	call checkNeighborList(nGroups, q_group, skinThickness, &
879		update_nlist)
880		#endif
881		endif
#	Line 750 | Line 926 | contains
926		endif
927
928		#ifdef IS_MPI
929	+	me_j = atid_col(j)
930		call get_interatomic_vector(q_group_Row(:,i), &
931		q_group_Col(:,j), d_grp, rgrpsq)
932		#else
933	+	me_j = atid(j)
934		call get_interatomic_vector(q_group(:,i), &
935		q_group(:,j), d_grp, rgrpsq)
936	<	#endif
936	>	#endif
937
938	<	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
938	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
939		if (update_nlist) then
940		nlist = nlist + 1
941
#	Line 777 | Line 955 | contains
955
956		list(nlist) = j
957		endif
958	+
959
960	<	if (loop .eq. PAIR_LOOP) then
961	<	vij = 0.0d0
783	<	fij(1:3) = 0.0d0
784	<	endif
960	>
961	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
962
963	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
964	<	in_switching_region)
965	<
966	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
967	<
968	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
969	<
970	<	atom1 = groupListRow(ia)
971	<
972	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
973	<
974	<	atom2 = groupListCol(jb)
975	<
976	<	if (skipThisPair(atom1, atom2)) cycle inner
977	<
978	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
979	<	d_atm(1:3) = d_grp(1:3)
980	<	ratmsq = rgrpsq
981	<	else
982	<	#ifdef IS_MPI
983	<	call get_interatomic_vector(q_Row(:,atom1), &
984	<	q_Col(:,atom2), d_atm, ratmsq)
985	<	#else
986	<	call get_interatomic_vector(q(:,atom1), &
987	<	q(:,atom2), d_atm, ratmsq)
988	<	#endif
989	<	endif
990	<
991	<	if (loop .eq. PREPAIR_LOOP) then
963	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
964	>	if (loop .eq. PAIR_LOOP) then
965	>	vij = 0.0d0
966	>	fij(1:3) = 0.0d0
967	>	endif
968	>
969	>	call get_switch(rgrpsq, sw, dswdr, rgrp, &
970	>	group_switch, in_switching_region)
971	>
972	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
973	>
974	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
975	>
976	>	atom1 = groupListRow(ia)
977	>
978	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
979	>
980	>	atom2 = groupListCol(jb)
981	>
982	>	if (skipThisPair(atom1, atom2))  cycle inner
983	>
984	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
985	>	d_atm(1:3) = d_grp(1:3)
986	>	ratmsq = rgrpsq
987	>	else
988	>	#ifdef IS_MPI
989	>	call get_interatomic_vector(q_Row(:,atom1), &
990	>	q_Col(:,atom2), d_atm, ratmsq)
991	>	#else
992	>	call get_interatomic_vector(q(:,atom1), &
993	>	q(:,atom2), d_atm, ratmsq)
994	>	#endif
995	>	endif
996	>
997	>	if (loop .eq. PREPAIR_LOOP) then
998		#ifdef IS_MPI
999	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1000	<	rgrpsq, d_grp, do_pot, do_stress, &
1001	<	eFrame, A, f, t, pot_local)
999	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1000	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1001	>	eFrame, A, f, t, pot_local)
1002		#else
1003	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1004	<	rgrpsq, d_grp, do_pot, do_stress, &
1005	<	eFrame, A, f, t, pot)
1003	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1004	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1005	>	eFrame, A, f, t, pot)
1006		#endif
1007	<	else
1007	>	else
1008		#ifdef IS_MPI
1009	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1010	<	do_pot, &
1011	<	eFrame, A, f, t, pot_local, vpair, fpair)
1009	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1010	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1011	>	fpair, d_grp, rgrp, rCut)
1012		#else
1013	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1014	<	do_pot,  &
1015	<	eFrame, A, f, t, pot, vpair, fpair)
1013	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1014	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1015	>	d_grp, rgrp, rCut)
1016		#endif
1017	+	vij = vij + vpair
1018	+	fij(1:3) = fij(1:3) + fpair(1:3)
1019	+	endif
1020	+	enddo inner
1021	+	enddo
1022
1023	<	vij = vij + vpair
1024	<	fij(1:3) = fij(1:3) + fpair(1:3)
1025	<	endif
1026	<	enddo inner
1027	<	enddo
1028	<
1029	<	if (loop .eq. PAIR_LOOP) then
1030	<	if (in_switching_region) then
1031	<	swderiv = vij*dswdr/rgrp
1032	<	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)
1023	>	if (loop .eq. PAIR_LOOP) then
1024	>	if (in_switching_region) then
1025	>	swderiv = vij*dswdr/rgrp
1026	>	fij(1) = fij(1) + swderiv*d_grp(1)
1027	>	fij(2) = fij(2) + swderiv*d_grp(2)
1028	>	fij(3) = fij(3) + swderiv*d_grp(3)
1029	>
1030	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1031	>	atom1=groupListRow(ia)
1032	>	mf = mfactRow(atom1)
1033		#ifdef IS_MPI
1034	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1035	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1036	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1034	>	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1035	>	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1036	>	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1037		#else
1038	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1039	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1040	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1038	>	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1039	>	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1040	>	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1041		#endif
1042	<	enddo
1043	<
1044	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1045	<	atom2=groupListCol(jb)
1046	<	mf = mfactCol(atom2)
1042	>	enddo
1043	>
1044	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1045	>	atom2=groupListCol(jb)
1046	>	mf = mfactCol(atom2)
1047		#ifdef IS_MPI
1048	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1049	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1050	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1048	>	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1049	>	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1050	>	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1051		#else
1052	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1053	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1054	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1052	>	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1053	>	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1054	>	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1055		#endif
1056	<	enddo
1057	<	endif
1056	>	enddo
1057	>	endif
1058
1059	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1059	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1060	>	endif
1061		endif
1062	<	end if
1062	>	endif
1063		enddo
1064	+
1065		enddo outer
1066
1067		if (update_nlist) then
#	Line 937 | Line 1121 | contains
1121
1122		if (do_pot) then
1123		! scatter/gather pot_row into the members of my column
1124	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1125	<
1124	>	do i = 1,LR_POT_TYPES
1125	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1126	>	end do
1127		! scatter/gather pot_local into all other procs
1128		! add resultant to get total pot
1129		do i = 1, nlocal
1130	<	pot_local = pot_local + pot_Temp(i)
1130	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1131	>	+ pot_Temp(1:LR_POT_TYPES,i)
1132		enddo
1133
1134		pot_Temp = 0.0_DP
1135	<
1136	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1135	>	do i = 1,LR_POT_TYPES
1136	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1137	>	end do
1138		do i = 1, nlocal
1139	<	pot_local = pot_local + pot_Temp(i)
1139	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1140	>	+ pot_Temp(1:LR_POT_TYPES,i)
1141		enddo
1142
1143		endif
1144		#endif
1145
1146	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1146	>	if (SIM_requires_postpair_calc) then
1147	>	do i = 1, nlocal
1148	>
1149	>	! we loop only over the local atoms, so we don't need row and column
1150	>	! lookups for the types
1151	>
1152	>	me_i = atid(i)
1153	>
1154	>	! is the atom electrostatic?  See if it would have an
1155	>	! electrostatic interaction with itself
1156	>	iHash = InteractionHash(me_i,me_i)
1157
1158	<	if (FF_uses_RF .and. SIM_uses_RF) then
961	<
1158	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1159		#ifdef IS_MPI
1160	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1161	<	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)
1160	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1161	>	t, do_pot)
1162		#else
1163	<	me_i = atid(i)
1163	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1164	>	t, do_pot)
1165		#endif
1166	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1166	>	endif
1167	>
1168	>
1169	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1170
1171	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1172	<
1173	<	mu_i = getDipoleMoment(me_i)
1174	<
1175	<	!! The reaction field needs to include a self contribution
1176	<	!! to the field:
1177	<	call accumulate_self_rf(i, mu_i, eFrame)
1178	<	!! Get the reaction field contribution to the
1179	<	!! potential and torques:
1180	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1171	>	! loop over the excludes to accumulate RF stuff we've
1172	>	! left out of the normal pair loop
1173	>
1174	>	do i1 = 1, nSkipsForAtom(i)
1175	>	j = skipsForAtom(i, i1)
1176	>
1177	>	! prevent overcounting of the skips
1178	>	if (i.lt.j) then
1179	>	call get_interatomic_vector(q(:,i), &
1180	>	q(:,j), d_atm, ratmsq)
1181	>	rVal = dsqrt(ratmsq)
1182	>	call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1183	>	in_switching_region)
1184		#ifdef IS_MPI
1185	<	pot_local = pot_local + rfpot
1185	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1186	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1187		#else
1188	<	pot = pot + rfpot
1189	<
1188	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1189	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1190		#endif
1191	<	endif
1192	<	enddo
1193	<	endif
1191	>	endif
1192	>	enddo
1193	>	endif
1194	>	enddo
1195		endif
1196	<
1001	<
1196	>
1197		#ifdef IS_MPI
1198	<
1198	>
1199		if (do_pot) then
1200	<	pot = pot + pot_local
1201	<	!! 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...
1200	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1201	>	mpi_comm_world,mpi_err)
1202		endif
1203	<
1203	>
1204		if (do_stress) then
1205		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1206		mpi_comm_world,mpi_err)
1207		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1208		mpi_comm_world,mpi_err)
1209		endif
1210	<
1210	>
1211		#else
1212	<
1212	>
1213		if (do_stress) then
1214		tau = tau_Temp
1215		virial = virial_Temp
1216		endif
1217	<
1217	>
1218		#endif
1219	<
1219	>
1220		end subroutine do_force_loop
1221
1222		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1223	<	eFrame, A, f, t, pot, vpair, fpair)
1223	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1224
1225	<	real( kind = dp ) :: pot, vpair, sw
1225	>	real( kind = dp ) :: vpair, sw
1226	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1227		real( kind = dp ), dimension(3) :: fpair
1228		real( kind = dp ), dimension(nLocal)   :: mfact
1229		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1039 | Line 1234 | contains
1234		logical, intent(inout) :: do_pot
1235		integer, intent(in) :: i, j
1236		real ( kind = dp ), intent(inout) :: rijsq
1237	<	real ( kind = dp )                :: r
1237	>	real ( kind = dp ), intent(inout) :: r_grp
1238		real ( kind = dp ), intent(inout) :: d(3)
1239	<	real ( kind = dp ) :: ebalance
1239	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1240	>	real ( kind = dp ), intent(inout) :: rCut
1241	>	real ( kind = dp ) :: r
1242		integer :: me_i, me_j
1243
1244	<	integer :: iMap
1244	>	integer :: iHash
1245
1246		r = sqrt(rijsq)
1247		vpair = 0.0d0
#	Line 1058 | Line 1255 | contains
1255		me_j = atid(j)
1256		#endif
1257
1258	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1259	<
1260	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1261	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1262	<	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)
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	<
1258	>	iHash = InteractionHash(me_i, me_j)
1259	>
1260	>	if ( iand(iHash, LJ_PAIR).ne.0 ) then
1261	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1262	>	pot(VDW_POT), f, do_pot)
1263		endif
1264	<
1265	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1264	>
1265	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1266	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1267	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1268	>	endif
1269	>
1270	>	if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1271		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1272	<	pot, A, f, t, do_pot)
1272	>	pot(HB_POT), A, f, t, do_pot)
1273		endif
1274	<
1275	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1274	>
1275	>	if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1276		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1277	<	pot, A, f, t, do_pot)
1277	>	pot(HB_POT), A, f, t, do_pot)
1278		endif
1279	<
1280	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1279	>
1280	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1281		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1282	<	pot, A, f, t, do_pot)
1282	>	pot(VDW_POT), A, f, t, do_pot)
1283		endif
1284
1285	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1286	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1287	<	!           pot, A, f, t, do_pot)
1285	>	if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1286	>	call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1287	>	pot(VDW_POT), A, f, t, do_pot)
1288		endif
1289	<
1290	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1291	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1292	<	do_pot)
1289	>
1290	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1291	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1292	>	pot(METALLIC_POT), f, do_pot)
1293		endif
1294	<
1295	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1294	>
1295	>	if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1296		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1297	<	pot, A, f, t, do_pot)
1297	>	pot(VDW_POT), A, f, t, do_pot)
1298		endif
1299	<
1300	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1299	>
1300	>	if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1301		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1302	<	pot, A, f, t, do_pot)
1302	>	pot(VDW_POT), A, f, t, do_pot)
1303		endif
1304	+
1305	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1306	+	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1307	+	pot(METALLIC_POT), f, do_pot)
1308	+	endif
1309	+
1310
1311	+
1312		end subroutine do_pair
1313
1314	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1314	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1315		do_pot, do_stress, eFrame, A, f, t, pot)
1316
1317	<	real( kind = dp ) :: pot, sw
1317	>	real( kind = dp ) :: sw
1318	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1319		real( kind = dp ), dimension(9,nLocal) :: eFrame
1320		real (kind=dp), dimension(9,nLocal) :: A
1321		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1125 | Line 1323 | contains
1323
1324		logical, intent(inout) :: do_pot, do_stress
1325		integer, intent(in) :: i, j
1326	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1326	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1327		real ( kind = dp )                :: r, rc
1328		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1329
1330	<	integer :: me_i, me_j, iMap
1330	>	integer :: me_i, me_j, iHash
1331
1332	+	r = sqrt(rijsq)
1333	+
1334		#ifdef IS_MPI
1335		me_i = atid_row(i)
1336		me_j = atid_col(j)
#	Line 1139 | Line 1339 | contains
1339		me_j = atid(j)
1340		#endif
1341
1342	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1342	>	iHash = InteractionHash(me_i, me_j)
1343
1344	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1345	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1344	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1345	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1346	>	endif
1347	>
1348	>	if ( iand(iHash, SC_PAIR).ne.0 ) then
1349	>	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1350		endif
1351
1352		end subroutine do_prepair
#	Line 1150 | Line 1354 | contains
1354
1355		subroutine do_preforce(nlocal,pot)
1356		integer :: nlocal
1357	<	real( kind = dp ) :: pot
1357	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1358
1359		if (FF_uses_EAM .and. SIM_uses_EAM) then
1360	<	call calc_EAM_preforce_Frho(nlocal,pot)
1360	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1361		endif
1362	+	if (FF_uses_SC .and. SIM_uses_SC) then
1363	+	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1364	+	endif
1365
1366
1367		end subroutine do_preforce
#	Line 1239 | Line 1446 | contains
1446		pot_Col = 0.0_dp
1447		pot_Temp = 0.0_dp
1448
1242	–	rf_Row = 0.0_dp
1243	–	rf_Col = 0.0_dp
1244	–	rf_Temp = 0.0_dp
1245	–
1449		#endif
1450
1451		if (FF_uses_EAM .and. SIM_uses_EAM) then
1452		call clean_EAM()
1453		endif
1454
1252	–	rf = 0.0_dp
1455		tau_Temp = 0.0_dp
1456		virial_Temp = 0.0_dp
1457		end subroutine zero_work_arrays
#	Line 1338 | Line 1540 | contains
1540
1541		function FF_UsesDirectionalAtoms() result(doesit)
1542		logical :: doesit
1543	<	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
1543	>	doesit = FF_uses_DirectionalAtoms
1544		end function FF_UsesDirectionalAtoms
1545
1546		function FF_RequiresPrepairCalc() result(doesit)
1547		logical :: doesit
1548	<	doesit = FF_uses_EAM
1548	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1549	>	.or. FF_uses_MEAM
1550		end function FF_RequiresPrepairCalc
1551
1351	–	function FF_RequiresPostpairCalc() result(doesit)
1352	–	logical :: doesit
1353	–	doesit = FF_uses_RF
1354	–	end function FF_RequiresPostpairCalc
1355	–
1552		#ifdef PROFILE
1553		function getforcetime() result(totalforcetime)
1554		real(kind=dp) :: totalforcetime

Diff Legend

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