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 2301 by gezelter, Thu Sep 15 22:05:21 2005 UTC vs.
Revision 2756 by gezelter, Wed May 17 15:37:15 2006 UTC

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

Diff Legend

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