ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE-4/src/UseTheForce/doForces.F90

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2267 by tim, Fri Jul 29 17:34:06 2005 UTC vs.
Revision 2533 by chuckv, Fri Dec 30 23:15:59 2005 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.25 2005-07-29 17:34:06 tim Exp $, $Date: 2005-07-29 17:34:06 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
48	>	!! @version $Id: doForces.F90,v 1.74 2005-12-30 23:15:59 chuckv Exp $, $Date: 2005-12-30 23:15:59 $, $Name: not supported by cvs2svn $, $Revision: 1.74 $
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, intent(out) :: status
156		integer :: i
157		integer :: j
158	<	integer :: ihash
157	<	real(kind=dp) :: myRcut
158	>	integer :: iHash
159		!! Test Types
160		logical :: i_is_LJ
161		logical :: i_is_Elect
#	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	<	status = 0
178	<
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!")
178	<	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 198 | 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 211 | 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 232 | 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 241 | 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
271
272	<	haveInteractionMap = .true.
273	<	end subroutine createInteractionMap
272	>	haveInteractionHash = .true.
273	>	end subroutine createInteractionHash
274
275	<	! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff.
255	<	subroutine createRcuts(defaultRList,stat)
256	<	real(kind=dp), intent(in), optional :: defaultRList
257	<	integer :: iMap
258	<	integer :: map_i,map_j
259	<	real(kind=dp) :: thisRCut = 0.0_dp
260	<	real(kind=dp) :: actualCutoff = 0.0_dp
261	<	integer, intent(out) :: stat
262	<	integer :: nAtypes
263	<	integer :: myStatus
275	>	subroutine createGtypeCutoffMap()
276
277	<	stat = 0
278	<	if (.not. haveInteractionMap) then
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	<	call createInteractionMap(myStatus)
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 (myStatus .ne. 0) then
296	<	write(default_error, *) 'createInteractionMap failed in doForces!'
297	<	stat = -1
298	<	return
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 (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
354	>	biggestAtypeCutoff = atypeMaxCutoff(i)
355	>	endif
356	>
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
392	<	nAtypes = getSize(atypes)
393	<	!! If we pass a default rcut, set all atypes to that cutoff distance
394	<	if(present(defaultRList)) then
395	<	InteractionMap(:,:)%rList = defaultRList
396	<	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
397	<	haveRlist = .true.
398	<	return
392	>
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	<	do map_i = 1,nAtypes
403	<	do map_j = map_i,nAtypes
404	<	iMap = InteractionMap(map_i, map_j)%InteractionHash
405	<
406	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
407	<	! thisRCut = getLJCutOff(map_i,map_j)
408	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
409	<	endif
410	<
411	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
412	<	! thisRCut = getElectrostaticCutOff(map_i,map_j)
413	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
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 ( iand(iMap, STICKY_PAIR).ne.0 ) then
302	<	! thisRCut = getStickyCutOff(map_i,map_j)
303	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
304	<	endif
305	<
306	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
307	<	! thisRCut = getStickyPowerCutOff(map_i,map_j)
308	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
309	<	endif
310	<
311	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
312	<	! thisRCut = getGayberneCutOff(map_i,map_j)
313	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
314	<	endif
315	<
316	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
317	<	!              thisRCut = getGaybrneLJCutOff(map_i,map_j)
318	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
319	<	endif
320	<
321	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
322	<	!              thisRCut = getEAMCutOff(map_i,map_j)
323	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
324	<	endif
325	<
326	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
327	<	!              thisRCut = getShapeCutOff(map_i,map_j)
328	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
329	<	endif
330	<
331	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
332	<	!              thisRCut = getShapeLJCutOff(map_i,map_j)
333	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
334	<	endif
335	<	InteractionMap(map_i, map_j)%rList = actualCutoff
336	<	InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
337	<	end do
338	<	end do
339	<	haveRlist = .true.
340	<	end subroutine createRcuts
528	>	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
529
530	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
531
532	<	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
533	<	!!$  subroutine setRlistDF( this_rlist )
534	<	!!$
346	<	!!$   real(kind=dp) :: this_rlist
347	<	!!$
348	<	!!$    rlist = this_rlist
349	<	!!$    rlistsq = rlist * rlist
350	<	!!$
351	<	!!$    haveRlist = .true.
352	<	!!$
353	<	!!$  end subroutine setRlistDF
532	>	if (.not.haveSkinThickness) then
533	>	skinThickness = 1.0_dp
534	>	endif
535
536	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
537
538	<	subroutine setSimVariables()
357	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
358	<	SIM_uses_LennardJones = SimUsesLennardJones()
359	<	SIM_uses_Electrostatics = SimUsesElectrostatics()
360	<	SIM_uses_Charges = SimUsesCharges()
361	<	SIM_uses_Dipoles = SimUsesDipoles()
362	<	SIM_uses_Sticky = SimUsesSticky()
363	<	SIM_uses_StickyPower = SimUsesStickyPower()
364	<	SIM_uses_GayBerne = SimUsesGayBerne()
365	<	SIM_uses_EAM = SimUsesEAM()
366	<	SIM_uses_Shapes = SimUsesShapes()
367	<	SIM_uses_FLARB = SimUsesFLARB()
368	<	SIM_uses_RF = SimUsesRF()
369	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
370	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
371	<	SIM_uses_PBC = SimUsesPBC()
538	>	! sanity check
539
540	<	haveSIMvariables = .true.
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	<	return
549	<	end subroutine setSimVariables
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	+
594	+	haveDefaultCutoffs = .true.
595	+	haveGtypeCutoffMap = .false.
596	+	end subroutine setCutoffs
597	+
598	+	subroutine cWasLame()
599	+
600	+	VisitCutoffsAfterComputing = .true.
601	+	return
602	+
603	+	end subroutine cWasLame
604	+
605	+	subroutine setCutoffPolicy(cutPolicy)
606	+
607	+	integer, intent(in) :: cutPolicy
608	+
609	+	cutoffPolicy = cutPolicy
610	+	haveCutoffPolicy = .true.
611	+	haveGtypeCutoffMap = .false.
612	+
613	+	end subroutine setCutoffPolicy
614	+
615	+	subroutine setElectrostaticMethod( thisESM )
616	+
617	+	integer, intent(in) :: thisESM
618	+
619	+	electrostaticSummationMethod = thisESM
620	+	haveElectrostaticSummationMethod = .true.
621	+
622	+	end subroutine setElectrostaticMethod
623	+
624	+	subroutine setSkinThickness( thisSkin )
625	+
626	+	real(kind=dp), intent(in) :: thisSkin
627	+
628	+	skinThickness = thisSkin
629	+	haveSkinThickness = .true.
630	+	haveGtypeCutoffMap = .false.
631	+
632	+	end subroutine setSkinThickness
633	+
634	+	subroutine setSimVariables()
635	+	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
636	+	SIM_uses_EAM = SimUsesEAM()
637	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
638	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
639	+	SIM_uses_PBC = SimUsesPBC()
640	+
641	+	haveSIMvariables = .true.
642	+
643	+	return
644	+	end subroutine setSimVariables
645	+
646		subroutine doReadyCheck(error)
647		integer, intent(out) :: error
648
#	Line 382 | Line 650 | contains
650
651		error = 0
652
653	<	if (.not. haveInteractionMap) then
653	>	if (.not. haveInteractionHash) then
654	>	call createInteractionHash()
655	>	endif
656
657	<	myStatus = 0
657	>	if (.not. haveGtypeCutoffMap) then
658	>	call createGtypeCutoffMap()
659	>	endif
660
389	–	call createInteractionMap(myStatus)
661
662	<	if (myStatus .ne. 0) then
663	<	write(default_error, *) 'createInteractionMap failed in doForces!'
393	<	error = -1
394	<	return
395	<	endif
662	>	if (VisitCutoffsAfterComputing) then
663	>	call set_switch(GROUP_SWITCH, largestRcut, largestRcut)
664		endif
665
666	+
667		if (.not. haveSIMvariables) then
668		call setSimVariables()
669		endif
670
671	<	if (.not. haveRlist) then
672	<	write(default_error, *) 'rList has not been set in doForces!'
673	<	error = -1
674	<	return
675	<	endif
671	>	!  if (.not. haveRlist) then
672	>	!     write(default_error, *) 'rList has not been set in doForces!'
673	>	!     error = -1
674	>	!     return
675	>	!  endif
676
677		if (.not. haveNeighborList) then
678		write(default_error, *) 'neighbor list has not been initialized in doForces!'
#	Line 428 | Line 697 | contains
697		end subroutine doReadyCheck
698
699
700	<	subroutine init_FF(use_RF_c, thisStat)
700	>	subroutine init_FF(thisStat)
701
433	–	logical, intent(in) :: use_RF_c
434	–
702		integer, intent(out) :: thisStat
703		integer :: my_status, nMatches
704		integer, pointer :: MatchList(:) => null()
438	–	real(kind=dp) :: rcut, rrf, rt, dielect
705
706		!! assume things are copacetic, unless they aren't
707		thisStat = 0
442	–
443	–	!! Fortran's version of a cast:
444	–	FF_uses_RF = use_RF_c
708
709		!! init_FF is called *after* all of the atom types have been
710		!! defined in atype_module using the new_atype subroutine.
#	Line 450 | Line 713 | contains
713		!! interactions are used by the force field.
714
715		FF_uses_DirectionalAtoms = .false.
453	–	FF_uses_LennardJones = .false.
454	–	FF_uses_Electrostatics = .false.
455	–	FF_uses_Charges = .false.
716		FF_uses_Dipoles = .false.
457	–	FF_uses_Sticky = .false.
458	–	FF_uses_StickyPower = .false.
717		FF_uses_GayBerne = .false.
718		FF_uses_EAM = .false.
719	<	FF_uses_Shapes = .false.
462	<	FF_uses_FLARB = .false.
719	>	FF_uses_SC = .false.
720
721		call getMatchingElementList(atypes, "is_Directional", .true., &
722		nMatches, MatchList)
723		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
724
468	–	call getMatchingElementList(atypes, "is_LennardJones", .true., &
469	–	nMatches, MatchList)
470	–	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
471	–
472	–	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
473	–	nMatches, MatchList)
474	–	if (nMatches .gt. 0) then
475	–	FF_uses_Electrostatics = .true.
476	–	endif
477	–
478	–	call getMatchingElementList(atypes, "is_Charge", .true., &
479	–	nMatches, MatchList)
480	–	if (nMatches .gt. 0) then
481	–	FF_uses_Charges = .true.
482	–	FF_uses_Electrostatics = .true.
483	–	endif
484	–
725		call getMatchingElementList(atypes, "is_Dipole", .true., &
726		nMatches, MatchList)
727	<	if (nMatches .gt. 0) then
488	<	FF_uses_Dipoles = .true.
489	<	FF_uses_Electrostatics = .true.
490	<	FF_uses_DirectionalAtoms = .true.
491	<	endif
492	<
493	<	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
494	<	nMatches, MatchList)
495	<	if (nMatches .gt. 0) then
496	<	FF_uses_Quadrupoles = .true.
497	<	FF_uses_Electrostatics = .true.
498	<	FF_uses_DirectionalAtoms = .true.
499	<	endif
500	<
501	<	call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
502	<	MatchList)
503	<	if (nMatches .gt. 0) then
504	<	FF_uses_Sticky = .true.
505	<	FF_uses_DirectionalAtoms = .true.
506	<	endif
507	<
508	<	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
509	<	MatchList)
510	<	if (nMatches .gt. 0) then
511	<	FF_uses_StickyPower = .true.
512	<	FF_uses_DirectionalAtoms = .true.
513	<	endif
727	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
728
729		call getMatchingElementList(atypes, "is_GayBerne", .true., &
730		nMatches, MatchList)
731	<	if (nMatches .gt. 0) then
518	<	FF_uses_GayBerne = .true.
519	<	FF_uses_DirectionalAtoms = .true.
520	<	endif
731	>	if (nMatches .gt. 0) FF_uses_GayBerne = .true.
732
733		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
734		if (nMatches .gt. 0) FF_uses_EAM = .true.
735
736	<	call getMatchingElementList(atypes, "is_Shape", .true., &
737	<	nMatches, MatchList)
527	<	if (nMatches .gt. 0) then
528	<	FF_uses_Shapes = .true.
529	<	FF_uses_DirectionalAtoms = .true.
530	<	endif
736	>	call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
737	>	if (nMatches .gt. 0) FF_uses_SC = .true.
738
532	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
533	–	nMatches, MatchList)
534	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
739
536	–	!! Assume sanity (for the sake of argument)
740		haveSaneForceField = .true.
538	–
539	–	!! check to make sure the FF_uses_RF setting makes sense
540	–
541	–	if (FF_uses_dipoles) then
542	–	if (FF_uses_RF) then
543	–	dielect = getDielect()
544	–	call initialize_rf(dielect)
545	–	endif
546	–	else
547	–	if (FF_uses_RF) then
548	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
549	–	thisStat = -1
550	–	haveSaneForceField = .false.
551	–	return
552	–	endif
553	–	endif
741
555	–	!sticky module does not contain check_sticky_FF anymore
556	–	!if (FF_uses_sticky) then
557	–	!   call check_sticky_FF(my_status)
558	–	!   if (my_status /= 0) then
559	–	!      thisStat = -1
560	–	!      haveSaneForceField = .false.
561	–	!      return
562	–	!   end if
563	–	!endif
564	–
742		if (FF_uses_EAM) then
743		call init_EAM_FF(my_status)
744		if (my_status /= 0) then
#	Line 572 | Line 749 | contains
749		end if
750		endif
751
575	–	if (FF_uses_GayBerne) then
576	–	call check_gb_pair_FF(my_status)
577	–	if (my_status .ne. 0) then
578	–	thisStat = -1
579	–	haveSaneForceField = .false.
580	–	return
581	–	endif
582	–	endif
583	–
584	–	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
585	–	endif
586	–
752		if (.not. haveNeighborList) then
753		!! Create neighbor lists
754		call expandNeighborList(nLocal, my_status)
#	Line 617 | Line 782 | contains
782
783		!! Stress Tensor
784		real( kind = dp), dimension(9) :: tau
785	<	real ( kind = dp ) :: pot
785	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
786		logical ( kind = 2) :: do_pot_c, do_stress_c
787		logical :: do_pot
788		logical :: do_stress
789		logical :: in_switching_region
790		#ifdef IS_MPI
791	<	real( kind = DP ) :: pot_local
791	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
792		integer :: nAtomsInRow
793		integer :: nAtomsInCol
794		integer :: nprocs
#	Line 638 | Line 803 | contains
803		integer :: nlist
804		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
805		real( kind = DP ) :: sw, dswdr, swderiv, mf
806	+	real( kind = DP ) :: rVal
807		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
808		real(kind=dp) :: rfpot, mu_i, virial
809	+	real(kind=dp):: rCut
810		integer :: me_i, me_j, n_in_i, n_in_j
811		logical :: is_dp_i
812		integer :: neighborListSize
#	Line 647 | Line 814 | contains
814		integer :: localError
815		integer :: propPack_i, propPack_j
816		integer :: loopStart, loopEnd, loop
817	<	integer :: iMap
818	<	real(kind=dp) :: listSkin = 1.0
817	>	integer :: iHash
818	>	integer :: i1
819	>
820
821		!! initialize local variables
822
#	Line 712 | Line 880 | contains
880		! (but only on the first time through):
881		if (loop .eq. loopStart) then
882		#ifdef IS_MPI
883	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
883	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
884		update_nlist)
885		#else
886	<	call checkNeighborList(nGroups, q_group, listSkin, &
886	>	call checkNeighborList(nGroups, q_group, skinThickness, &
887		update_nlist)
888		#endif
889		endif
#	Line 739 | Line 907 | contains
907		#endif
908		outer: do i = istart, iend
909
742	–	#ifdef IS_MPI
743	–	me_i = atid_row(i)
744	–	#else
745	–	me_i = atid(i)
746	–	#endif
747	–
910		if (update_nlist) point(i) = nlist + 1
911
912		n_in_i = groupStartRow(i+1) - groupStartRow(i)
#	Line 779 | Line 941 | contains
941		me_j = atid(j)
942		call get_interatomic_vector(q_group(:,i), &
943		q_group(:,j), d_grp, rgrpsq)
944	<	#endif
944	>	#endif
945
946	<	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
946	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
947		if (update_nlist) then
948		nlist = nlist + 1
949
#	Line 801 | Line 963 | contains
963
964		list(nlist) = j
965		endif
966	+
967
968	<	if (loop .eq. PAIR_LOOP) then
969	<	vij = 0.0d0
807	<	fij(1:3) = 0.0d0
808	<	endif
968	>
969	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
970
971	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
972	<	in_switching_region)
973	<
974	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
975	<
976	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
977	<
978	<	atom1 = groupListRow(ia)
979	<
980	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
981	<
982	<	atom2 = groupListCol(jb)
983	<
984	<	if (skipThisPair(atom1, atom2)) cycle inner
985	<
986	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
987	<	d_atm(1:3) = d_grp(1:3)
988	<	ratmsq = rgrpsq
989	<	else
990	<	#ifdef IS_MPI
991	<	call get_interatomic_vector(q_Row(:,atom1), &
992	<	q_Col(:,atom2), d_atm, ratmsq)
971	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
972	>	if (loop .eq. PAIR_LOOP) then
973	>	vij = 0.0d0
974	>	fij(1:3) = 0.0d0
975	>	endif
976	>
977	>	call get_switch(rgrpsq, sw, dswdr, rgrp, &
978	>	group_switch, in_switching_region)
979	>
980	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
981	>
982	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
983	>
984	>	atom1 = groupListRow(ia)
985	>
986	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
987	>
988	>	atom2 = groupListCol(jb)
989	>
990	>	if (skipThisPair(atom1, atom2))  cycle inner
991	>
992	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
993	>	d_atm(1:3) = d_grp(1:3)
994	>	ratmsq = rgrpsq
995	>	else
996	>	#ifdef IS_MPI
997	>	call get_interatomic_vector(q_Row(:,atom1), &
998	>	q_Col(:,atom2), d_atm, ratmsq)
999		#else
1000	<	call get_interatomic_vector(q(:,atom1), &
1001	<	q(:,atom2), d_atm, ratmsq)
1000	>	call get_interatomic_vector(q(:,atom1), &
1001	>	q(:,atom2), d_atm, ratmsq)
1002		#endif
1003	<	endif
1004	<
1005	<	if (loop .eq. PREPAIR_LOOP) then
1003	>	endif
1004	>
1005	>	if (loop .eq. PREPAIR_LOOP) then
1006		#ifdef IS_MPI
1007	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1008	<	rgrpsq, d_grp, do_pot, do_stress, &
1009	<	eFrame, A, f, t, pot_local)
1007	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1008	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1009	>	eFrame, A, f, t, pot_local)
1010		#else
1011	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1012	<	rgrpsq, d_grp, do_pot, do_stress, &
1013	<	eFrame, A, f, t, pot)
1011	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1012	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1013	>	eFrame, A, f, t, pot)
1014		#endif
1015	<	else
1015	>	else
1016		#ifdef IS_MPI
1017	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1018	<	do_pot, &
1019	<	eFrame, A, f, t, pot_local, vpair, fpair)
1017	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1018	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1019	>	fpair, d_grp, rgrp, rCut)
1020		#else
1021	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1022	<	do_pot,  &
1023	<	eFrame, A, f, t, pot, vpair, fpair)
1021	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1022	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1023	>	d_grp, rgrp, rCut)
1024		#endif
1025	+	vij = vij + vpair
1026	+	fij(1:3) = fij(1:3) + fpair(1:3)
1027	+	endif
1028	+	enddo inner
1029	+	enddo
1030
1031	<	vij = vij + vpair
1032	<	fij(1:3) = fij(1:3) + fpair(1:3)
1033	<	endif
1034	<	enddo inner
1035	<	enddo
1036	<
1037	<	if (loop .eq. PAIR_LOOP) then
1038	<	if (in_switching_region) then
1039	<	swderiv = vij*dswdr/rgrp
1040	<	fij(1) = fij(1) + swderiv*d_grp(1)
869	<	fij(2) = fij(2) + swderiv*d_grp(2)
870	<	fij(3) = fij(3) + swderiv*d_grp(3)
871	<
872	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
873	<	atom1=groupListRow(ia)
874	<	mf = mfactRow(atom1)
1031	>	if (loop .eq. PAIR_LOOP) then
1032	>	if (in_switching_region) then
1033	>	swderiv = vij*dswdr/rgrp
1034	>	fij(1) = fij(1) + swderiv*d_grp(1)
1035	>	fij(2) = fij(2) + swderiv*d_grp(2)
1036	>	fij(3) = fij(3) + swderiv*d_grp(3)
1037	>
1038	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1039	>	atom1=groupListRow(ia)
1040	>	mf = mfactRow(atom1)
1041		#ifdef IS_MPI
1042	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1043	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1044	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1042	>	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1043	>	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1044	>	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1045		#else
1046	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1047	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1048	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1046	>	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1047	>	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1048	>	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1049		#endif
1050	<	enddo
1051	<
1052	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1053	<	atom2=groupListCol(jb)
1054	<	mf = mfactCol(atom2)
1050	>	enddo
1051	>
1052	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1053	>	atom2=groupListCol(jb)
1054	>	mf = mfactCol(atom2)
1055		#ifdef IS_MPI
1056	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1057	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1058	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1056	>	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1057	>	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1058	>	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1059		#else
1060	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1061	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1062	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1060	>	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1061	>	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1062	>	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1063		#endif
1064	<	enddo
1065	<	endif
1064	>	enddo
1065	>	endif
1066
1067	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1067	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1068	>	endif
1069		endif
1070	<	end if
1070	>	endif
1071		enddo
1072	+
1073		enddo outer
1074
1075		if (update_nlist) then
#	Line 961 | Line 1129 | contains
1129
1130		if (do_pot) then
1131		! scatter/gather pot_row into the members of my column
1132	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1133	<
1132	>	do i = 1,LR_POT_TYPES
1133	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1134	>	end do
1135		! scatter/gather pot_local into all other procs
1136		! add resultant to get total pot
1137		do i = 1, nlocal
1138	<	pot_local = pot_local + pot_Temp(i)
1138	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1139	>	+ pot_Temp(1:LR_POT_TYPES,i)
1140		enddo
1141
1142		pot_Temp = 0.0_DP
1143	<
1144	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1143	>	do i = 1,LR_POT_TYPES
1144	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1145	>	end do
1146		do i = 1, nlocal
1147	<	pot_local = pot_local + pot_Temp(i)
1147	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1148	>	+ pot_Temp(1:LR_POT_TYPES,i)
1149		enddo
1150
1151		endif
1152		#endif
1153
1154	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1154	>	if (SIM_requires_postpair_calc) then
1155	>	do i = 1, nlocal
1156	>
1157	>	! we loop only over the local atoms, so we don't need row and column
1158	>	! lookups for the types
1159	>
1160	>	me_i = atid(i)
1161	>
1162	>	! is the atom electrostatic?  See if it would have an
1163	>	! electrostatic interaction with itself
1164	>	iHash = InteractionHash(me_i,me_i)
1165
1166	<	if (FF_uses_RF .and. SIM_uses_RF) then
985	<
1166	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1167		#ifdef IS_MPI
1168	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1169	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
989	<	do i = 1,nlocal
990	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
991	<	end do
992	<	#endif
993	<
994	<	do i = 1, nLocal
995	<
996	<	rfpot = 0.0_DP
997	<	#ifdef IS_MPI
998	<	me_i = atid_row(i)
1168	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1169	>	t, do_pot)
1170		#else
1171	<	me_i = atid(i)
1171	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1172	>	t, do_pot)
1173		#endif
1174	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1174	>	endif
1175	>
1176	>
1177	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1178
1179	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1180	<
1181	<	mu_i = getDipoleMoment(me_i)
1182	<
1183	<	!! The reaction field needs to include a self contribution
1184	<	!! to the field:
1185	<	call accumulate_self_rf(i, mu_i, eFrame)
1186	<	!! Get the reaction field contribution to the
1187	<	!! potential and torques:
1188	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1179	>	! loop over the excludes to accumulate RF stuff we've
1180	>	! left out of the normal pair loop
1181	>
1182	>	do i1 = 1, nSkipsForAtom(i)
1183	>	j = skipsForAtom(i, i1)
1184	>
1185	>	! prevent overcounting of the skips
1186	>	if (i.lt.j) then
1187	>	call get_interatomic_vector(q(:,i), &
1188	>	q(:,j), d_atm, ratmsq)
1189	>	rVal = dsqrt(ratmsq)
1190	>	call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1191	>	in_switching_region)
1192		#ifdef IS_MPI
1193	<	pot_local = pot_local + rfpot
1193	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1194	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1195		#else
1196	<	pot = pot + rfpot
1197	<
1196	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1197	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1198		#endif
1199	<	endif
1200	<	enddo
1201	<	endif
1199	>	endif
1200	>	enddo
1201	>	endif
1202	>	enddo
1203		endif
1204	<
1025	<
1204	>
1205		#ifdef IS_MPI
1206	<
1206	>
1207		if (do_pot) then
1208	<	pot = pot + pot_local
1209	<	!! we assume the c code will do the allreduce to get the total potential
1031	<	!! we could do it right here if we needed to...
1208	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1209	>	mpi_comm_world,mpi_err)
1210		endif
1211	<
1211	>
1212		if (do_stress) then
1213		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1214		mpi_comm_world,mpi_err)
1215		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1216		mpi_comm_world,mpi_err)
1217		endif
1218	<
1218	>
1219		#else
1220	<
1220	>
1221		if (do_stress) then
1222		tau = tau_Temp
1223		virial = virial_Temp
1224		endif
1225	<
1225	>
1226		#endif
1227	<
1227	>
1228		end subroutine do_force_loop
1229
1230		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1231	<	eFrame, A, f, t, pot, vpair, fpair)
1231	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1232
1233	<	real( kind = dp ) :: pot, vpair, sw
1233	>	real( kind = dp ) :: vpair, sw
1234	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1235		real( kind = dp ), dimension(3) :: fpair
1236		real( kind = dp ), dimension(nLocal)   :: mfact
1237		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1063 | Line 1242 | contains
1242		logical, intent(inout) :: do_pot
1243		integer, intent(in) :: i, j
1244		real ( kind = dp ), intent(inout) :: rijsq
1245	<	real ( kind = dp )                :: r
1245	>	real ( kind = dp ), intent(inout) :: r_grp
1246		real ( kind = dp ), intent(inout) :: d(3)
1247	<	real ( kind = dp ) :: ebalance
1247	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1248	>	real ( kind = dp ), intent(inout) :: rCut
1249	>	real ( kind = dp ) :: r
1250		integer :: me_i, me_j
1251
1252	<	integer :: iMap
1252	>	integer :: iHash
1253
1254		r = sqrt(rijsq)
1255		vpair = 0.0d0
#	Line 1082 | Line 1263 | contains
1263		me_j = atid(j)
1264		#endif
1265
1266	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1267	<
1268	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1269	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1270	<	endif
1090	<
1091	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1092	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1093	<	pot, eFrame, f, t, do_pot)
1094	<
1095	<	if (FF_uses_RF .and. SIM_uses_RF) then
1096	<
1097	<	! CHECK ME (RF needs to know about all electrostatic types)
1098	<	call accumulate_rf(i, j, r, eFrame, sw)
1099	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1100	<	endif
1101	<
1266	>	iHash = InteractionHash(me_i, me_j)
1267	>
1268	>	if ( iand(iHash, LJ_PAIR).ne.0 ) then
1269	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1270	>	pot(VDW_POT), f, do_pot)
1271		endif
1272	<
1273	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1272	>
1273	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1274	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1275	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1276	>	endif
1277	>
1278	>	if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1279		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1280	<	pot, A, f, t, do_pot)
1280	>	pot(HB_POT), A, f, t, do_pot)
1281		endif
1282	<
1283	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1282	>
1283	>	if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1284		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1285	<	pot, A, f, t, do_pot)
1285	>	pot(HB_POT), A, f, t, do_pot)
1286		endif
1287	<
1288	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1287	>
1288	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1289		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1290	<	pot, A, f, t, do_pot)
1290	>	pot(VDW_POT), A, f, t, do_pot)
1291		endif
1292
1293	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1294	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1295	<	!           pot, A, f, t, do_pot)
1293	>	if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1294	>	call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1295	>	pot(VDW_POT), A, f, t, do_pot)
1296		endif
1297	<
1298	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1299	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1300	<	do_pot)
1297	>
1298	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1299	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1300	>	pot(METALLIC_POT), f, do_pot)
1301		endif
1302	<
1303	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1302	>
1303	>	if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1304		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1305	<	pot, A, f, t, do_pot)
1305	>	pot(VDW_POT), A, f, t, do_pot)
1306		endif
1307	<
1308	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1307	>
1308	>	if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1309		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1310	<	pot, A, f, t, do_pot)
1310	>	pot(VDW_POT), A, f, t, do_pot)
1311		endif
1312	+
1313	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1314	+	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1315	+	pot(METALLIC_POT), f, do_pot)
1316	+	endif
1317	+
1318
1319	+
1320		end subroutine do_pair
1321
1322	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1322	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1323		do_pot, do_stress, eFrame, A, f, t, pot)
1324
1325	<	real( kind = dp ) :: pot, sw
1325	>	real( kind = dp ) :: sw
1326	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1327		real( kind = dp ), dimension(9,nLocal) :: eFrame
1328		real (kind=dp), dimension(9,nLocal) :: A
1329		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1149 | Line 1331 | contains
1331
1332		logical, intent(inout) :: do_pot, do_stress
1333		integer, intent(in) :: i, j
1334	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1334	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1335		real ( kind = dp )                :: r, rc
1336		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1337
1338	<	integer :: me_i, me_j, iMap
1338	>	integer :: me_i, me_j, iHash
1339
1340	+	r = sqrt(rijsq)
1341	+
1342		#ifdef IS_MPI
1343		me_i = atid_row(i)
1344		me_j = atid_col(j)
#	Line 1163 | Line 1347 | contains
1347		me_j = atid(j)
1348		#endif
1349
1350	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1350	>	iHash = InteractionHash(me_i, me_j)
1351
1352	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1353	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1352	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1353	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1354		endif
1355	+
1356	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1357	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1358	+	endif
1359
1360		end subroutine do_prepair
1361
1362
1363		subroutine do_preforce(nlocal,pot)
1364		integer :: nlocal
1365	<	real( kind = dp ) :: pot
1365	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1366
1367		if (FF_uses_EAM .and. SIM_uses_EAM) then
1368	<	call calc_EAM_preforce_Frho(nlocal,pot)
1368	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1369		endif
1370	+	if (FF_uses_SC .and. SIM_uses_SC) then
1371	+	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1372	+	endif
1373
1374
1375		end subroutine do_preforce
#	Line 1263 | Line 1454 | contains
1454		pot_Col = 0.0_dp
1455		pot_Temp = 0.0_dp
1456
1266	–	rf_Row = 0.0_dp
1267	–	rf_Col = 0.0_dp
1268	–	rf_Temp = 0.0_dp
1269	–
1457		#endif
1458
1459		if (FF_uses_EAM .and. SIM_uses_EAM) then
1460		call clean_EAM()
1461		endif
1462
1276	–	rf = 0.0_dp
1463		tau_Temp = 0.0_dp
1464		virial_Temp = 0.0_dp
1465		end subroutine zero_work_arrays
#	Line 1362 | Line 1548 | contains
1548
1549		function FF_UsesDirectionalAtoms() result(doesit)
1550		logical :: doesit
1551	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1366	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1367	<	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1551	>	doesit = FF_uses_DirectionalAtoms
1552		end function FF_UsesDirectionalAtoms
1553
1554		function FF_RequiresPrepairCalc() result(doesit)
1555		logical :: doesit
1556	<	doesit = FF_uses_EAM
1556	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1557	>	.or. FF_uses_MEAM
1558		end function FF_RequiresPrepairCalc
1559
1375	–	function FF_RequiresPostpairCalc() result(doesit)
1376	–	logical :: doesit
1377	–	doesit = FF_uses_RF
1378	–	end function FF_RequiresPostpairCalc
1379	–
1560		#ifdef PROFILE
1561		function getforcetime() result(totalforcetime)
1562		real(kind=dp) :: totalforcetime

Diff Legend

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