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 2309 by chrisfen, Sun Sep 18 20:45:38 2005 UTC vs.
Revision 2758 by gezelter, Wed May 17 19:54:27 2006 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.46 2005-09-18 20:45:38 chrisfen Exp $, $Date: 2005-09-18 20:45:38 $, $Name: not supported by cvs2svn $, $Revision: 1.46 $
48	>	!! @version $Id: doForces.F90,v 1.82 2006-05-17 19:54:26 gezelter Exp $, $Date: 2006-05-17 19:54:26 $, $Name: not supported by cvs2svn $, $Revision: 1.82 $
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_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_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 :: 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
106	–	public :: createInteractionHash
107	–	public :: createGtypeCutoffMap
108	–	public :: getStickyCut
109	–	public :: getStickyPowerCut
110	–	public :: getGayBerneCut
111	–	public :: getEAMCut
112	–	public :: getShapeCut
123
124		#ifdef PROFILE
125		public :: getforcetime
#	Line 122 | 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 132 | Line 147 | module doForces
147		end type gtypeCutoffs
148		type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
149
150	<	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
136	<	real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
137	<
150	>
151		contains
152
153	<	subroutine createInteractionHash(status)
153	>	subroutine createInteractionHash()
154		integer :: nAtypes
142	–	integer, intent(out) :: status
155		integer :: i
156		integer :: j
157		integer :: iHash
#	Line 151 | 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 158 | 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
163	–	status = 0
164	–
179		if (.not. associated(atypes)) then
180	<	call handleError("atype", "atypes was not present before call of createInteractionHash!")
167	<	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
#	Line 191 | 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 204 | 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 225 | 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 244 | Line 271 | contains
271		haveInteractionHash = .true.
272		end subroutine createInteractionHash
273
274	<	subroutine createGtypeCutoffMap(stat)
274	>	subroutine createGtypeCutoffMap()
275
249	–	integer, intent(out), optional :: stat
276		logical :: i_is_LJ
277		logical :: i_is_Elect
278		logical :: i_is_Sticky
#	Line 254 | 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
265	–	stat = 0
294		if (.not. haveInteractionHash) then
295	<	call createInteractionHash(myStatus)
268	<	if (myStatus .ne. 0) then
269	<	write(default_error, *) 'createInteractionHash failed in doForces!'
270	<	stat = -1
271	<	return
272	<	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 286 | 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
325	–
326	–	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 356 | 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
405	–	skin = defaultRlist - defaultRcut
406	–	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)
415	<	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
416	<	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	<	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()
609	<	end subroutine setCutoffPolicy
430	<
608	>	haveCutoffPolicy = .true.
609	>	haveGtypeCutoffMap = .false.
610
611	<	subroutine setSimVariables()
612	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
613	<	SIM_uses_EAM = SimUsesEAM()
435	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
436	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
437	<	SIM_uses_PBC = SimUsesPBC()
611	>	end subroutine setCutoffPolicy
612	>
613	>	subroutine setElectrostaticMethod( thisESM )
614
615	<	haveSIMvariables = .true.
615	>	integer, intent(in) :: thisESM
616
617	<	return
618	<	end subroutine setSimVariables
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	+	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	+
645		subroutine doReadyCheck(error)
646		integer, intent(out) :: error
446	–
647		integer :: myStatus
648
649		error = 0
650
651		if (.not. haveInteractionHash) then
652	<	myStatus = 0
453	<	call createInteractionHash(myStatus)
454	<	if (myStatus .ne. 0) then
455	<	write(default_error, *) 'createInteractionHash failed in doForces!'
456	<	error = -1
457	<	return
458	<	endif
652	>	call createInteractionHash()
653		endif
654
655		if (.not. haveGtypeCutoffMap) then
656	<	myStatus = 0
463	<	call createGtypeCutoffMap(myStatus)
464	<	if (myStatus .ne. 0) then
465	<	write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
466	<	error = -1
467	<	return
468	<	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
475	–	!  if (.not. haveRlist) then
476	–	!     write(default_error, *) 'rList has not been set in doForces!'
477	–	!     error = -1
478	–	!     return
479	–	!  endif
480	–
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 501 | Line 691 | contains
691		end subroutine doReadyCheck
692
693
694	<	subroutine init_FF(thisESM, thisStat)
694	>	subroutine init_FF(thisStat)
695
506	–	integer, intent(in) :: thisESM
507	–	real(kind=dp), intent(in) :: dampingAlpha
696		integer, intent(out) :: thisStat
697		integer :: my_status, nMatches
698		integer, pointer :: MatchList(:) => null()
511	–	real(kind=dp) :: rcut, rrf, rt, dielect
699
700		!! assume things are copacetic, unless they aren't
701		thisStat = 0
702
516	–	electrostaticSummationMethod = thisESM
517	–
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 525 | 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 541 | 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
547	–	!! check to make sure the reaction field setting makes sense
548	–
549	–	if (FF_uses_Dipoles) then
550	–	if (electrostaticSummationMethod == 3) then
551	–	dielect = getDielect()
552	–	call initialize_rf(dielect)
553	–	endif
554	–	else
555	–	if (electrostaticSummationMethod == 3) then
556	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
557	–	thisStat = -1
558	–	haveSaneForceField = .false.
559	–	return
560	–	endif
561	–	endif
562	–
736		if (FF_uses_EAM) then
737		call init_EAM_FF(my_status)
738		if (my_status /= 0) then
#	Line 570 | Line 743 | contains
743		end if
744		endif
745
573	–	if (FF_uses_GayBerne) then
574	–	call check_gb_pair_FF(my_status)
575	–	if (my_status .ne. 0) then
576	–	thisStat = -1
577	–	haveSaneForceField = .false.
578	–	return
579	–	endif
580	–	endif
581	–
746		if (.not. haveNeighborList) then
747		!! Create neighbor lists
748		call expandNeighborList(nLocal, my_status)
#	Line 612 | 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 633 | 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 643 | 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 707 | 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 768 | 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 790 | Line 957 | contains
957
958		list(nlist) = j
959		endif
960	+
961	+	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
962
963	<	if (loop .eq. PAIR_LOOP) then
964	<	vij = 0.0d0
965	<	fij(1:3) = 0.0d0
966	<	endif
967	<
968	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
969	<	in_switching_region)
970	<
971	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
972	<
973	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
974	<
975	<	atom1 = groupListRow(ia)
976	<
977	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
978	<
979	<	atom2 = groupListCol(jb)
980	<
981	<	if (skipThisPair(atom1, atom2)) cycle inner
982	<
983	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
984	<	d_atm(1:3) = d_grp(1:3)
985	<	ratmsq = rgrpsq
986	<	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)
858	<	fij(2) = fij(2) + swderiv*d_grp(2)
859	<	fij(3) = fij(3) + swderiv*d_grp(3)
860	<
861	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
862	<	atom1=groupListRow(ia)
863	<	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 950 | 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 (electrostaticSummationMethod == 3) then
974	<
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)
978	<	do i = 1,nlocal
979	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
980	<	end do
981	<	#endif
982	<
983	<	do i = 1, nLocal
984	<
985	<	rfpot = 0.0_DP
986	<	#ifdef IS_MPI
987	<	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	<
1014	<
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
1205	<	!! we could do it right here if we needed to...
1203	>	#ifdef SINGLE_PRECISION
1204	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1205	>	mpi_comm_world,mpi_err)
1206	>	#else
1207	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1208	>	mpi_comm_world,mpi_err)
1209	>	#endif
1210		endif
1211	<
1211	>
1212		if (do_stress) then
1213	+	#ifdef SINGLE_PRECISION
1214	+	call mpi_allreduce(tau_Temp, tau, 9,mpi_real,mpi_sum, &
1215	+	mpi_comm_world,mpi_err)
1216	+	call mpi_allreduce(virial_Temp, virial,1,mpi_real,mpi_sum, &
1217	+	mpi_comm_world,mpi_err)
1218	+	#else
1219		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1220		mpi_comm_world,mpi_err)
1221		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1222		mpi_comm_world,mpi_err)
1223	+	#endif
1224		endif
1225	<
1225	>
1226		#else
1227	<
1227	>
1228		if (do_stress) then
1229		tau = tau_Temp
1230		virial = virial_Temp
1231		endif
1232	<
1232	>
1233		#endif
1234	<
1234	>
1235		end subroutine do_force_loop
1236
1237		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1238	<	eFrame, A, f, t, pot, vpair, fpair)
1238	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1239
1240	<	real( kind = dp ) :: pot, vpair, sw
1240	>	real( kind = dp ) :: vpair, sw
1241	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1242		real( kind = dp ), dimension(3) :: fpair
1243		real( kind = dp ), dimension(nLocal)   :: mfact
1244		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1052 | Line 1249 | contains
1249		logical, intent(inout) :: do_pot
1250		integer, intent(in) :: i, j
1251		real ( kind = dp ), intent(inout) :: rijsq
1252	<	real ( kind = dp )                :: r
1252	>	real ( kind = dp ), intent(inout) :: r_grp
1253		real ( kind = dp ), intent(inout) :: d(3)
1254	+	real ( kind = dp ), intent(inout) :: d_grp(3)
1255	+	real ( kind = dp ), intent(inout) :: rCut
1256	+	real ( kind = dp ) :: r
1257	+	real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1258		integer :: me_i, me_j
1259	+	integer :: k
1260
1261		integer :: iHash
1262
1263		r = sqrt(rijsq)
1264	<	vpair = 0.0d0
1265	<	fpair(1:3) = 0.0d0
1264	>
1265	>	vpair = 0.0_dp
1266	>	fpair(1:3) = 0.0_dp
1267
1268		#ifdef IS_MPI
1269		me_i = atid_row(i)
#	Line 1071 | Line 1274 | contains
1274		#endif
1275
1276		iHash = InteractionHash(me_i, me_j)
1277	<
1277	>
1278		if ( iand(iHash, LJ_PAIR).ne.0 ) then
1279	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1279	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1280	>	pot(VDW_POT), f, do_pot)
1281		endif
1282	<
1282	>
1283		if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1284	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1285	<	pot, eFrame, f, t, do_pot)
1082	<
1083	<	if (electrostaticSummationMethod == 3) then
1084	<
1085	<	! CHECK ME (RF needs to know about all electrostatic types)
1086	<	call accumulate_rf(i, j, r, eFrame, sw)
1087	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1088	<	endif
1089	<
1284	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1285	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1286		endif
1287	<
1287	>
1288		if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1289		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1290	<	pot, A, f, t, do_pot)
1290	>	pot(HB_POT), A, f, t, do_pot)
1291		endif
1292	<
1292	>
1293		if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1294		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1295	<	pot, A, f, t, do_pot)
1295	>	pot(HB_POT), A, f, t, do_pot)
1296		endif
1297	<
1297	>
1298		if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1299		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1300	<	pot, A, f, t, do_pot)
1300	>	pot(VDW_POT), A, f, t, do_pot)
1301		endif
1302
1303		if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1304	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1305	<	!           pot, A, f, t, do_pot)
1304	>	call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1305	>	pot(VDW_POT), A, f, t, do_pot)
1306		endif
1307	<
1307	>
1308		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1309	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1310	<	do_pot)
1309	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1310	>	pot(METALLIC_POT), f, do_pot)
1311		endif
1312	<
1312	>
1313		if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1314		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1315	<	pot, A, f, t, do_pot)
1315	>	pot(VDW_POT), A, f, t, do_pot)
1316		endif
1317	<
1317	>
1318		if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1319		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1320	<	pot, A, f, t, do_pot)
1320	>	pot(VDW_POT), A, f, t, do_pot)
1321		endif
1322	<
1322	>
1323	>	if ( iand(iHash, SC_PAIR).ne.0 ) then
1324	>	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1325	>	pot(METALLIC_POT), f, do_pot)
1326	>	endif
1327	>
1328		end subroutine do_pair
1329
1330	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1330	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1331		do_pot, do_stress, eFrame, A, f, t, pot)
1332
1333	<	real( kind = dp ) :: pot, sw
1333	>	real( kind = dp ) :: sw
1334	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1335		real( kind = dp ), dimension(9,nLocal) :: eFrame
1336		real (kind=dp), dimension(9,nLocal) :: A
1337		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1137 | Line 1339 | contains
1339
1340		logical, intent(inout) :: do_pot, do_stress
1341		integer, intent(in) :: i, j
1342	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1342	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1343		real ( kind = dp )                :: r, rc
1344		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1345
1346		integer :: me_i, me_j, iHash
1347
1348		r = sqrt(rijsq)
1349	<
1349	>
1350		#ifdef IS_MPI
1351		me_i = atid_row(i)
1352		me_j = atid_col(j)
#	Line 1156 | Line 1358 | contains
1358		iHash = InteractionHash(me_i, me_j)
1359
1360		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1361	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1361	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1362		endif
1363	+
1364	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1365	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1366	+	endif
1367
1368		end subroutine do_prepair
1369
1370
1371		subroutine do_preforce(nlocal,pot)
1372		integer :: nlocal
1373	<	real( kind = dp ) :: pot
1373	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1374
1375		if (FF_uses_EAM .and. SIM_uses_EAM) then
1376	<	call calc_EAM_preforce_Frho(nlocal,pot)
1376	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1377		endif
1378	<
1379	<
1378	>	if (FF_uses_SC .and. SIM_uses_SC) then
1379	>	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1380	>	endif
1381		end subroutine do_preforce
1382
1383
#	Line 1182 | Line 1389 | contains
1389		real( kind = dp ) :: d(3), scaled(3)
1390		integer i
1391
1392	<	d(1:3) = q_j(1:3) - q_i(1:3)
1392	>	d(1) = q_j(1) - q_i(1)
1393	>	d(2) = q_j(2) - q_i(2)
1394	>	d(3) = q_j(3) - q_i(3)
1395
1396		! Wrap back into periodic box if necessary
1397		if ( SIM_uses_PBC ) then
1398
1399		if( .not.boxIsOrthorhombic ) then
1400		! calc the scaled coordinates.
1401	+	! scaled = matmul(HmatInv, d)
1402
1403	<	scaled = matmul(HmatInv, d)
1404	<
1403	>	scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1404	>	scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1405	>	scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1406	>
1407		! wrap the scaled coordinates
1408
1409	<	scaled = scaled  - anint(scaled)
1409	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1410	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1411	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1412
1199	–
1413		! calc the wrapped real coordinates from the wrapped scaled
1414		! coordinates
1415	<
1416	<	d = matmul(Hmat,scaled)
1415	>	! d = matmul(Hmat,scaled)
1416	>	d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1417	>	d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1418	>	d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1419
1420		else
1421		! calc the scaled coordinates.
1422
1423	<	do i = 1, 3
1424	<	scaled(i) = d(i) * HmatInv(i,i)
1423	>	scaled(1) = d(1) * HmatInv(1,1)
1424	>	scaled(2) = d(2) * HmatInv(2,2)
1425	>	scaled(3) = d(3) * HmatInv(3,3)
1426	>
1427	>	! wrap the scaled coordinates
1428	>
1429	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1430	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1431	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1432
1433	<	! wrap the scaled coordinates
1433	>	! calc the wrapped real coordinates from the wrapped scaled
1434	>	! coordinates
1435
1436	<	scaled(i) = scaled(i) - anint(scaled(i))
1436	>	d(1) = scaled(1)*Hmat(1,1)
1437	>	d(2) = scaled(2)*Hmat(2,2)
1438	>	d(3) = scaled(3)*Hmat(3,3)
1439
1215	–	! calc the wrapped real coordinates from the wrapped scaled
1216	–	! coordinates
1217	–
1218	–	d(i) = scaled(i)*Hmat(i,i)
1219	–	enddo
1440		endif
1441
1442		endif
1443
1444	<	r_sq = dot_product(d,d)
1444	>	r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1445
1446		end subroutine get_interatomic_vector
1447
#	Line 1253 | Line 1473 | contains
1473		pot_Col = 0.0_dp
1474		pot_Temp = 0.0_dp
1475
1256	–	rf_Row = 0.0_dp
1257	–	rf_Col = 0.0_dp
1258	–	rf_Temp = 0.0_dp
1259	–
1476		#endif
1477
1478		if (FF_uses_EAM .and. SIM_uses_EAM) then
1479		call clean_EAM()
1480		endif
1481
1266	–	rf = 0.0_dp
1482		tau_Temp = 0.0_dp
1483		virial_Temp = 0.0_dp
1484		end subroutine zero_work_arrays
#	Line 1357 | Line 1572 | contains
1572
1573		function FF_RequiresPrepairCalc() result(doesit)
1574		logical :: doesit
1575	<	doesit = FF_uses_EAM
1575	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1576	>	.or. FF_uses_MEAM
1577		end function FF_RequiresPrepairCalc
1578
1363	–	function FF_RequiresPostpairCalc() result(doesit)
1364	–	logical :: doesit
1365	–	if (electrostaticSummationMethod == 3) doesit = .true.
1366	–	end function FF_RequiresPostpairCalc
1367	–
1579		#ifdef PROFILE
1580		function getforcetime() result(totalforcetime)
1581		real(kind=dp) :: totalforcetime

Diff Legend

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