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 2722 by gezelter, Thu Apr 20 18:24:24 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.79 2006-04-20 18:24:24 gezelter Exp $, $Date: 2006-04-20 18:24:24 $, $Name: not supported by cvs2svn $, $Revision: 1.79 $
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	+	use interpolation
68		#ifdef IS_MPI
69		use mpiSimulation
70		#endif
#	Line 72 | Line 73 | module doForces
73		PRIVATE
74
75		#define __FORTRAN90
75	–	#include "UseTheForce/fSwitchingFunction.h"
76		#include "UseTheForce/fCutoffPolicy.h"
77		#include "UseTheForce/DarkSide/fInteractionMap.h"
78	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79
79	–
80		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81		INTEGER, PARAMETER:: PAIR_LOOP    = 2
82	+	INTEGER, PARAMETER:: np = 500
83
84		logical, save :: haveNeighborList = .false.
85		logical, save :: haveSIMvariables = .false.
86		logical, save :: haveSaneForceField = .false.
87		logical, save :: haveInteractionHash = .false.
88		logical, save :: haveGtypeCutoffMap = .false.
89	<	logical, save :: haveRlist = .false.
89	>	logical, save :: haveDefaultCutoffs = .false.
90	>	logical, save :: haveSkinThickness = .false.
91	>	logical, save :: haveElectrostaticSummationMethod = .false.
92	>	logical, save :: haveCutoffPolicy = .false.
93	>	logical, save :: VisitCutoffsAfterComputing = .false.
94	>	logical, save :: haveSplineSqrt = .false.
95
96		logical, save :: FF_uses_DirectionalAtoms
97		logical, save :: FF_uses_Dipoles
98		logical, save :: FF_uses_GayBerne
99		logical, save :: FF_uses_EAM
100	+	logical, save :: FF_uses_SC
101	+	logical, save :: FF_uses_MEAM
102	+
103
104		logical, save :: SIM_uses_DirectionalAtoms
105		logical, save :: SIM_uses_EAM
106	+	logical, save :: SIM_uses_SC
107	+	logical, save :: SIM_uses_MEAM
108		logical, save :: SIM_requires_postpair_calc
109		logical, save :: SIM_requires_prepair_calc
110		logical, save :: SIM_uses_PBC
111
112		integer, save :: electrostaticSummationMethod
113	+	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
114
115	+	real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
116	+	real(kind=dp), save :: skinThickness
117	+	logical, save :: defaultDoShift
118	+
119		public :: init_FF
120	<	public :: setDefaultCutoffs
120	>	public :: setCutoffs
121	>	public :: cWasLame
122	>	public :: setElectrostaticMethod
123	>	public :: setCutoffPolicy
124	>	public :: setSkinThickness
125		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
126
127		#ifdef PROFILE
128		public :: getforcetime
#	Line 122 | Line 135 | module doForces
135		! Bit hash to determine pair-pair interactions.
136		integer, dimension(:,:), allocatable :: InteractionHash
137		real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
138	<	real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
139	<	integer, dimension(:), allocatable :: groupToGtype
140	<	real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
138	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
139	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
140	>
141	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
142	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
143	>
144	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
145	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
146		type ::gtypeCutoffs
147		real(kind=dp) :: rcut
148		real(kind=dp) :: rcutsq
#	Line 132 | Line 150 | module doForces
150		end type gtypeCutoffs
151		type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
152
153	<	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
154	<	real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
153	>	! variables for the spline of the sqrt
154	>	type(cubicSpline), save :: splineSqrt
155	>	logical, save :: useSpline = .true.
156
157	+
158		contains
159
160	<	subroutine createInteractionHash(status)
160	>	subroutine createInteractionHash()
161		integer :: nAtypes
142	–	integer, intent(out) :: status
162		integer :: i
163		integer :: j
164		integer :: iHash
#	Line 151 | Line 170 | contains
170		logical :: i_is_GB
171		logical :: i_is_EAM
172		logical :: i_is_Shape
173	+	logical :: i_is_SC
174	+	logical :: i_is_MEAM
175		logical :: j_is_LJ
176		logical :: j_is_Elect
177		logical :: j_is_Sticky
#	Line 158 | Line 179 | contains
179		logical :: j_is_GB
180		logical :: j_is_EAM
181		logical :: j_is_Shape
182	+	logical :: j_is_SC
183	+	logical :: j_is_MEAM
184		real(kind=dp) :: myRcut
185
163	–	status = 0
164	–
186		if (.not. associated(atypes)) then
187	<	call handleError("atype", "atypes was not present before call of createInteractionHash!")
167	<	status = -1
187	>	call handleError("doForces", "atypes was not present before call of createInteractionHash!")
188		return
189		endif
190
191		nAtypes = getSize(atypes)
192
193		if (nAtypes == 0) then
194	<	status = -1
194	>	call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
195		return
196		end if
197
198		if (.not. allocated(InteractionHash)) then
199	+	allocate(InteractionHash(nAtypes,nAtypes))
200	+	else
201	+	deallocate(InteractionHash)
202		allocate(InteractionHash(nAtypes,nAtypes))
203		endif
204
205		if (.not. allocated(atypeMaxCutoff)) then
206		allocate(atypeMaxCutoff(nAtypes))
207	+	else
208	+	deallocate(atypeMaxCutoff)
209	+	allocate(atypeMaxCutoff(nAtypes))
210		endif
211
212		do i = 1, nAtypes
#	Line 191 | Line 217 | contains
217		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
218		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
219		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
220	+	call getElementProperty(atypes, i, "is_SC", i_is_SC)
221	+	call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
222
223		do j = i, nAtypes
224
#	Line 204 | Line 232 | contains
232		call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
233		call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
234		call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
235	+	call getElementProperty(atypes, j, "is_SC", j_is_SC)
236	+	call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
237
238		if (i_is_LJ .and. j_is_LJ) then
239		iHash = ior(iHash, LJ_PAIR)
#	Line 225 | Line 255 | contains
255		iHash = ior(iHash, EAM_PAIR)
256		endif
257
258	+	if (i_is_SC .and. j_is_SC) then
259	+	iHash = ior(iHash, SC_PAIR)
260	+	endif
261	+
262		if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
263		if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
264		if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
#	Line 244 | Line 278 | contains
278		haveInteractionHash = .true.
279		end subroutine createInteractionHash
280
281	<	subroutine createGtypeCutoffMap(stat)
281	>	subroutine createGtypeCutoffMap()
282
249	–	integer, intent(out), optional :: stat
283		logical :: i_is_LJ
284		logical :: i_is_Elect
285		logical :: i_is_Sticky
#	Line 254 | Line 287 | contains
287		logical :: i_is_GB
288		logical :: i_is_EAM
289		logical :: i_is_Shape
290	+	logical :: i_is_SC
291		logical :: GtypeFound
292
293		integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
294	<	integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
294	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
295		integer :: nGroupsInRow
296	<	real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
296	>	integer :: nGroupsInCol
297	>	integer :: nGroupTypesRow,nGroupTypesCol
298	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
299		real(kind=dp) :: biggestAtypeCutoff
300
265	–	stat = 0
301		if (.not. haveInteractionHash) then
302	<	call createInteractionHash(myStatus)
268	<	if (myStatus .ne. 0) then
269	<	write(default_error, *) 'createInteractionHash failed in doForces!'
270	<	stat = -1
271	<	return
272	<	endif
302	>	call createInteractionHash()
303		endif
304		#ifdef IS_MPI
305		nGroupsInRow = getNgroupsInRow(plan_group_row)
306	+	nGroupsInCol = getNgroupsInCol(plan_group_col)
307		#endif
308		nAtypes = getSize(atypes)
309		! Set all of the initial cutoffs to zero.
#	Line 286 | Line 317 | contains
317		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
318		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
319		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
320	<
320	>	call getElementProperty(atypes, i, "is_SC", i_is_SC)
321
322	<	if (i_is_LJ) then
323	<	thisRcut = getSigma(i) * 2.5_dp
324	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325	<	endif
326	<	if (i_is_Elect) then
327	<	thisRcut = defaultRcut
328	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329	<	endif
330	<	if (i_is_Sticky) then
331	<	thisRcut = getStickyCut(i)
332	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333	<	endif
334	<	if (i_is_StickyP) then
335	<	thisRcut = getStickyPowerCut(i)
336	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337	<	endif
338	<	if (i_is_GB) then
339	<	thisRcut = getGayBerneCut(i)
340	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
341	<	endif
342	<	if (i_is_EAM) then
343	<	thisRcut = getEAMCut(i)
344	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
345	<	endif
346	<	if (i_is_Shape) then
347	<	thisRcut = getShapeCut(i)
348	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322	>	if (haveDefaultCutoffs) then
323	>	atypeMaxCutoff(i) = defaultRcut
324	>	else
325	>	if (i_is_LJ) then
326	>	thisRcut = getSigma(i) * 2.5_dp
327	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
328	>	endif
329	>	if (i_is_Elect) then
330	>	thisRcut = defaultRcut
331	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
332	>	endif
333	>	if (i_is_Sticky) then
334	>	thisRcut = getStickyCut(i)
335	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
336	>	endif
337	>	if (i_is_StickyP) then
338	>	thisRcut = getStickyPowerCut(i)
339	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
340	>	endif
341	>	if (i_is_GB) then
342	>	thisRcut = getGayBerneCut(i)
343	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
344	>	endif
345	>	if (i_is_EAM) then
346	>	thisRcut = getEAMCut(i)
347	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
348	>	endif
349	>	if (i_is_Shape) then
350	>	thisRcut = getShapeCut(i)
351	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
352	>	endif
353	>	if (i_is_SC) then
354	>	thisRcut = getSCCut(i)
355	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
356	>	endif
357		endif
358	<
358	>
359		if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
360		biggestAtypeCutoff = atypeMaxCutoff(i)
361		endif
362	+
363		endif
364		enddo
325	–
326	–	nGroupTypes = 0
365
366		istart = 1
367	+	jstart = 1
368		#ifdef IS_MPI
369		iend = nGroupsInRow
370	+	jend = nGroupsInCol
371		#else
372		iend = nGroups
373	+	jend = nGroups
374		#endif
375
376		!! allocate the groupToGtype and gtypeMaxCutoff here.
377	<	if(.not.allocated(groupToGtype)) then
378	<	allocate(groupToGtype(iend))
379	<	allocate(groupMaxCutoff(iend))
380	<	allocate(gtypeMaxCutoff(iend))
381	<	groupMaxCutoff = 0.0_dp
382	<	gtypeMaxCutoff = 0.0_dp
377	>	if(.not.allocated(groupToGtypeRow)) then
378	>	!  allocate(groupToGtype(iend))
379	>	allocate(groupToGtypeRow(iend))
380	>	else
381	>	deallocate(groupToGtypeRow)
382	>	allocate(groupToGtypeRow(iend))
383		endif
384	+	if(.not.allocated(groupMaxCutoffRow)) then
385	+	allocate(groupMaxCutoffRow(iend))
386	+	else
387	+	deallocate(groupMaxCutoffRow)
388	+	allocate(groupMaxCutoffRow(iend))
389	+	end if
390	+
391	+	if(.not.allocated(gtypeMaxCutoffRow)) then
392	+	allocate(gtypeMaxCutoffRow(iend))
393	+	else
394	+	deallocate(gtypeMaxCutoffRow)
395	+	allocate(gtypeMaxCutoffRow(iend))
396	+	endif
397	+
398	+
399	+	#ifdef IS_MPI
400	+	! We only allocate new storage if we are in MPI because Ncol /= Nrow
401	+	if(.not.associated(groupToGtypeCol)) then
402	+	allocate(groupToGtypeCol(jend))
403	+	else
404	+	deallocate(groupToGtypeCol)
405	+	allocate(groupToGtypeCol(jend))
406	+	end if
407	+
408	+	if(.not.associated(groupMaxCutoffCol)) then
409	+	allocate(groupMaxCutoffCol(jend))
410	+	else
411	+	deallocate(groupMaxCutoffCol)
412	+	allocate(groupMaxCutoffCol(jend))
413	+	end if
414	+	if(.not.associated(gtypeMaxCutoffCol)) then
415	+	allocate(gtypeMaxCutoffCol(jend))
416	+	else
417	+	deallocate(gtypeMaxCutoffCol)
418	+	allocate(gtypeMaxCutoffCol(jend))
419	+	end if
420	+
421	+	groupMaxCutoffCol = 0.0_dp
422	+	gtypeMaxCutoffCol = 0.0_dp
423	+
424	+	#endif
425	+	groupMaxCutoffRow = 0.0_dp
426	+	gtypeMaxCutoffRow = 0.0_dp
427	+
428	+
429		!! first we do a single loop over the cutoff groups to find the
430		!! largest cutoff for any atypes present in this group.  We also
431		!! create gtypes at this point.
432
433		tol = 1.0d-6
434	<
434	>	nGroupTypesRow = 0
435	>	nGroupTypesCol = 0
436		do i = istart, iend
437		n_in_i = groupStartRow(i+1) - groupStartRow(i)
438	<	groupMaxCutoff(i) = 0.0_dp
438	>	groupMaxCutoffRow(i) = 0.0_dp
439		do ia = groupStartRow(i), groupStartRow(i+1)-1
440		atom1 = groupListRow(ia)
441		#ifdef IS_MPI
#	Line 356 | Line 443 | contains
443		#else
444		me_i = atid(atom1)
445		#endif
446	<	if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
447	<	groupMaxCutoff(i)=atypeMaxCutoff(me_i)
446	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
447	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
448		endif
449		enddo
450	+	if (nGroupTypesRow.eq.0) then
451	+	nGroupTypesRow = nGroupTypesRow + 1
452	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
453	+	groupToGtypeRow(i) = nGroupTypesRow
454	+	else
455	+	GtypeFound = .false.
456	+	do g = 1, nGroupTypesRow
457	+	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
458	+	groupToGtypeRow(i) = g
459	+	GtypeFound = .true.
460	+	endif
461	+	enddo
462	+	if (.not.GtypeFound) then
463	+	nGroupTypesRow = nGroupTypesRow + 1
464	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
465	+	groupToGtypeRow(i) = nGroupTypesRow
466	+	endif
467	+	endif
468	+	enddo
469
470	<	if (nGroupTypes.eq.0) then
471	<	nGroupTypes = nGroupTypes + 1
472	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
473	<	groupToGtype(i) = nGroupTypes
470	>	#ifdef IS_MPI
471	>	do j = jstart, jend
472	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
473	>	groupMaxCutoffCol(j) = 0.0_dp
474	>	do ja = groupStartCol(j), groupStartCol(j+1)-1
475	>	atom1 = groupListCol(ja)
476	>
477	>	me_j = atid_col(atom1)
478	>
479	>	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
480	>	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
481	>	endif
482	>	enddo
483	>
484	>	if (nGroupTypesCol.eq.0) then
485	>	nGroupTypesCol = nGroupTypesCol + 1
486	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
487	>	groupToGtypeCol(j) = nGroupTypesCol
488		else
489		GtypeFound = .false.
490	<	do g = 1, nGroupTypes
491	<	if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
492	<	groupToGtype(i) = g
490	>	do g = 1, nGroupTypesCol
491	>	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
492	>	groupToGtypeCol(j) = g
493		GtypeFound = .true.
494		endif
495		enddo
496		if (.not.GtypeFound) then
497	<	nGroupTypes = nGroupTypes + 1
498	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
499	<	groupToGtype(i) = nGroupTypes
497	>	nGroupTypesCol = nGroupTypesCol + 1
498	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
499	>	groupToGtypeCol(j) = nGroupTypesCol
500		endif
501		endif
502		enddo
503
504	+	#else
505	+	! Set pointers to information we just found
506	+	nGroupTypesCol = nGroupTypesRow
507	+	groupToGtypeCol => groupToGtypeRow
508	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
509	+	groupMaxCutoffCol => groupMaxCutoffRow
510	+	#endif
511	+
512		!! allocate the gtypeCutoffMap here.
513	<	allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
513	>	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
514		!! then we do a double loop over all the group TYPES to find the cutoff
515		!! map between groups of two types
516	<
517	<	do i = 1, nGroupTypes
518	<	do j = 1, nGroupTypes
516	>	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
517	>
518	>	do i = 1, nGroupTypesRow
519	>	do j = 1, nGroupTypesCol
520
521		select case(cutoffPolicy)
522		case(TRADITIONAL_CUTOFF_POLICY)
523	<	thisRcut = maxval(gtypeMaxCutoff)
523	>	thisRcut = tradRcut
524		case(MIX_CUTOFF_POLICY)
525	<	thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
525	>	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
526		case(MAX_CUTOFF_POLICY)
527	<	thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
527	>	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
528		case default
529		call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
530		return
531		end select
532		gtypeCutoffMap(i,j)%rcut = thisRcut
533	+
534	+	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
535	+
536		gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
405	–	skin = defaultRlist - defaultRcut
406	–	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
537
538	+	if (.not.haveSkinThickness) then
539	+	skinThickness = 1.0_dp
540	+	endif
541	+
542	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
543	+
544	+	! sanity check
545	+
546	+	if (haveDefaultCutoffs) then
547	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
548	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
549	+	endif
550	+	endif
551		enddo
552		enddo
553	+
554	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
555	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
556	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
557	+	#ifdef IS_MPI
558	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
559	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
560	+	#endif
561	+	groupMaxCutoffCol => null()
562	+	gtypeMaxCutoffCol => null()
563
564		haveGtypeCutoffMap = .true.
565		end subroutine createGtypeCutoffMap
566
567	<	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
415	<	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
416	<	integer, intent(in) :: cutPolicy
567	>	subroutine setCutoffs(defRcut, defRsw)
568
569	+	real(kind=dp),intent(in) :: defRcut, defRsw
570	+	character(len = statusMsgSize) :: errMsg
571	+	integer :: localError
572	+
573		defaultRcut = defRcut
574		defaultRsw = defRsw
575	<	defaultRlist = defRlist
576	<	cutoffPolicy = cutPolicy
577	<	end subroutine setDefaultCutoffs
575	>
576	>	defaultDoShift = .false.
577	>	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
578	>
579	>	write(errMsg, *) &
580	>	'cutoffRadius and switchingRadius are set to the same', newline &
581	>	// tab, 'value.  OOPSE will use shifted ', newline &
582	>	// tab, 'potentials instead of switching functions.'
583	>
584	>	call handleInfo("setCutoffs", errMsg)
585	>
586	>	defaultDoShift = .true.
587	>
588	>	endif
589	>
590	>	localError = 0
591	>	call setLJDefaultCutoff( defaultRcut, defaultDoShift )
592	>	call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
593	>	call setCutoffEAM( defaultRcut )
594	>	call setCutoffSC( defaultRcut )
595	>	call set_switch(defaultRsw, defaultRcut)
596	>	call setHmatDangerousRcutValue(defaultRcut)
597	>	call setupSplineSqrt(defaultRcut)
598	>
599	>	haveDefaultCutoffs = .true.
600	>	haveGtypeCutoffMap = .false.
601
602	<	subroutine setCutoffPolicy(cutPolicy)
602	>	end subroutine setCutoffs
603
604	+	subroutine cWasLame()
605	+
606	+	VisitCutoffsAfterComputing = .true.
607	+	return
608	+
609	+	end subroutine cWasLame
610	+
611	+	subroutine setCutoffPolicy(cutPolicy)
612	+
613		integer, intent(in) :: cutPolicy
614	+
615		cutoffPolicy = cutPolicy
616	<	call createGtypeCutoffMap()
617	<	end subroutine setCutoffPolicy
430	<
616	>	haveCutoffPolicy = .true.
617	>	haveGtypeCutoffMap = .false.
618
619	<	subroutine setSimVariables()
620	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
621	<	SIM_uses_EAM = SimUsesEAM()
435	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
436	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
437	<	SIM_uses_PBC = SimUsesPBC()
619	>	end subroutine setCutoffPolicy
620	>
621	>	subroutine setElectrostaticMethod( thisESM )
622
623	<	haveSIMvariables = .true.
623	>	integer, intent(in) :: thisESM
624
625	<	return
626	<	end subroutine setSimVariables
625	>	electrostaticSummationMethod = thisESM
626	>	haveElectrostaticSummationMethod = .true.
627	>
628	>	end subroutine setElectrostaticMethod
629	>
630	>	subroutine setSkinThickness( thisSkin )
631	>
632	>	real(kind=dp), intent(in) :: thisSkin
633	>
634	>	skinThickness = thisSkin
635	>	haveSkinThickness = .true.
636	>	haveGtypeCutoffMap = .false.
637	>
638	>	end subroutine setSkinThickness
639	>
640	>	subroutine setSimVariables()
641	>	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
642	>	SIM_uses_EAM = SimUsesEAM()
643	>	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
644	>	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
645	>	SIM_uses_PBC = SimUsesPBC()
646	>	SIM_uses_SC = SimUsesSC()
647	>
648	>	haveSIMvariables = .true.
649	>
650	>	return
651	>	end subroutine setSimVariables
652
653		subroutine doReadyCheck(error)
654		integer, intent(out) :: error
446	–
655		integer :: myStatus
656
657		error = 0
658
659		if (.not. haveInteractionHash) then
660	<	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
660	>	call createInteractionHash()
661		endif
662
663		if (.not. haveGtypeCutoffMap) then
664	<	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
664	>	call createGtypeCutoffMap()
665		endif
666
667	+	if (VisitCutoffsAfterComputing) then
668	+	call set_switch(largestRcut, largestRcut)
669	+	call setHmatDangerousRcutValue(largestRcut)
670	+	call setLJsplineRmax(largestRcut)
671	+	call setCutoffEAM(largestRcut)
672	+	call setCutoffSC(largestRcut)
673	+	VisitCutoffsAfterComputing = .false.
674	+	endif
675	+
676	+	if (.not. haveSplineSqrt) then
677	+	call setupSplineSqrt(largestRcut)
678	+	endif
679	+
680		if (.not. haveSIMvariables) then
681		call setSimVariables()
682		endif
683
475	–	!  if (.not. haveRlist) then
476	–	!     write(default_error, *) 'rList has not been set in doForces!'
477	–	!     error = -1
478	–	!     return
479	–	!  endif
480	–
684		if (.not. haveNeighborList) then
685		write(default_error, *) 'neighbor list has not been initialized in doForces!'
686		error = -1
687		return
688		end if
689	<
689	>
690		if (.not. haveSaneForceField) then
691		write(default_error, *) 'Force Field is not sane in doForces!'
692		error = -1
693		return
694		end if
695	<
695	>
696		#ifdef IS_MPI
697		if (.not. isMPISimSet()) then
698		write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
#	Line 501 | Line 704 | contains
704		end subroutine doReadyCheck
705
706
707	<	subroutine init_FF(thisESM, thisStat)
707	>	subroutine init_FF(thisStat)
708
506	–	integer, intent(in) :: thisESM
507	–	real(kind=dp), intent(in) :: dampingAlpha
709		integer, intent(out) :: thisStat
710		integer :: my_status, nMatches
711		integer, pointer :: MatchList(:) => null()
511	–	real(kind=dp) :: rcut, rrf, rt, dielect
712
713		!! assume things are copacetic, unless they aren't
714		thisStat = 0
715
516	–	electrostaticSummationMethod = thisESM
517	–
716		!! init_FF is called *after* all of the atom types have been
717		!! defined in atype_module using the new_atype subroutine.
718		!!
#	Line 525 | Line 723 | contains
723		FF_uses_Dipoles = .false.
724		FF_uses_GayBerne = .false.
725		FF_uses_EAM = .false.
726	+	FF_uses_SC = .false.
727
728		call getMatchingElementList(atypes, "is_Directional", .true., &
729		nMatches, MatchList)
#	Line 541 | Line 740 | contains
740		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
741		if (nMatches .gt. 0) FF_uses_EAM = .true.
742
743	+	call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
744	+	if (nMatches .gt. 0) FF_uses_SC = .true.
745
545	–	haveSaneForceField = .true.
746
747	<	!! check to make sure the reaction field setting makes sense
747	>	haveSaneForceField = .true.
748
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	–
749		if (FF_uses_EAM) then
750		call init_EAM_FF(my_status)
751		if (my_status /= 0) then
#	Line 570 | Line 756 | contains
756		end if
757		endif
758
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	–
759		if (.not. haveNeighborList) then
760		!! Create neighbor lists
761		call expandNeighborList(nLocal, my_status)
#	Line 612 | Line 789 | contains
789
790		!! Stress Tensor
791		real( kind = dp), dimension(9) :: tau
792	<	real ( kind = dp ) :: pot
792	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
793		logical ( kind = 2) :: do_pot_c, do_stress_c
794		logical :: do_pot
795		logical :: do_stress
796		logical :: in_switching_region
797		#ifdef IS_MPI
798	<	real( kind = DP ) :: pot_local
798	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
799		integer :: nAtomsInRow
800		integer :: nAtomsInCol
801		integer :: nprocs
#	Line 633 | Line 810 | contains
810		integer :: nlist
811		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
812		real( kind = DP ) :: sw, dswdr, swderiv, mf
813	+	real( kind = DP ) :: rVal
814		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
815		real(kind=dp) :: rfpot, mu_i, virial
816	+	real(kind=dp):: rCut
817		integer :: me_i, me_j, n_in_i, n_in_j
818		logical :: is_dp_i
819		integer :: neighborListSize
#	Line 643 | Line 822 | contains
822		integer :: propPack_i, propPack_j
823		integer :: loopStart, loopEnd, loop
824		integer :: iHash
825	<	real(kind=dp) :: listSkin = 1.0
825	>	integer :: i1
826	>
827
828		!! initialize local variables
829
#	Line 707 | Line 887 | contains
887		! (but only on the first time through):
888		if (loop .eq. loopStart) then
889		#ifdef IS_MPI
890	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
890	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
891		update_nlist)
892		#else
893	<	call checkNeighborList(nGroups, q_group, listSkin, &
893	>	call checkNeighborList(nGroups, q_group, skinThickness, &
894		update_nlist)
895		#endif
896		endif
#	Line 768 | Line 948 | contains
948		me_j = atid(j)
949		call get_interatomic_vector(q_group(:,i), &
950		q_group(:,j), d_grp, rgrpsq)
951	<	#endif
951	>	#endif
952
953	<	if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
953	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
954		if (update_nlist) then
955		nlist = nlist + 1
956
#	Line 790 | Line 970 | contains
970
971		list(nlist) = j
972		endif
973	+
974	+	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
975
976	<	if (loop .eq. PAIR_LOOP) then
977	<	vij = 0.0d0
978	<	fij(1:3) = 0.0d0
979	<	endif
980	<
981	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
982	<	in_switching_region)
983	<
984	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
985	<
986	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
987	<
988	<	atom1 = groupListRow(ia)
989	<
990	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
991	<
992	<	atom2 = groupListCol(jb)
993	<
994	<	if (skipThisPair(atom1, atom2)) cycle inner
995	<
996	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
997	<	d_atm(1:3) = d_grp(1:3)
998	<	ratmsq = rgrpsq
999	<	else
976	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
977	>	if (loop .eq. PAIR_LOOP) then
978	>	vij = 0.0d0
979	>	fij(1) = 0.0_dp
980	>	fij(2) = 0.0_dp
981	>	fij(3) = 0.0_dp
982	>	endif
983	>
984	>	call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
985	>
986	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
987	>
988	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
989	>
990	>	atom1 = groupListRow(ia)
991	>
992	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
993	>
994	>	atom2 = groupListCol(jb)
995	>
996	>	if (skipThisPair(atom1, atom2))  cycle inner
997	>
998	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
999	>	d_atm(1) = d_grp(1)
1000	>	d_atm(2) = d_grp(2)
1001	>	d_atm(3) = d_grp(3)
1002	>	ratmsq = rgrpsq
1003	>	else
1004		#ifdef IS_MPI
1005	<	call get_interatomic_vector(q_Row(:,atom1), &
1006	<	q_Col(:,atom2), d_atm, ratmsq)
1005	>	call get_interatomic_vector(q_Row(:,atom1), &
1006	>	q_Col(:,atom2), d_atm, ratmsq)
1007		#else
1008	<	call get_interatomic_vector(q(:,atom1), &
1009	<	q(:,atom2), d_atm, ratmsq)
1008	>	call get_interatomic_vector(q(:,atom1), &
1009	>	q(:,atom2), d_atm, ratmsq)
1010		#endif
1011	<	endif
1012	<
1013	<	if (loop .eq. PREPAIR_LOOP) then
1011	>	endif
1012	>
1013	>	if (loop .eq. PREPAIR_LOOP) then
1014		#ifdef IS_MPI
1015	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1016	<	rgrpsq, d_grp, do_pot, do_stress, &
1017	<	eFrame, A, f, t, pot_local)
1015	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1016	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1017	>	eFrame, A, f, t, pot_local)
1018		#else
1019	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1020	<	rgrpsq, d_grp, do_pot, do_stress, &
1021	<	eFrame, A, f, t, pot)
1019	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1020	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1021	>	eFrame, A, f, t, pot)
1022		#endif
1023	<	else
1023	>	else
1024		#ifdef IS_MPI
1025	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1026	<	do_pot, &
1027	<	eFrame, A, f, t, pot_local, vpair, fpair)
1025	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1026	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1027	>	fpair, d_grp, rgrp, rCut)
1028		#else
1029	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1030	<	do_pot,  &
1031	<	eFrame, A, f, t, pot, vpair, fpair)
1029	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1030	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1031	>	d_grp, rgrp, rCut)
1032		#endif
1033	+	vij = vij + vpair
1034	+	fij(1) = fij(1) + fpair(1)
1035	+	fij(2) = fij(2) + fpair(2)
1036	+	fij(3) = fij(3) + fpair(3)
1037	+	endif
1038	+	enddo inner
1039	+	enddo
1040
1041	<	vij = vij + vpair
1042	<	fij(1:3) = fij(1:3) + fpair(1:3)
1043	<	endif
1044	<	enddo inner
1045	<	enddo
1046	<
1047	<	if (loop .eq. PAIR_LOOP) then
1048	<	if (in_switching_region) then
1049	<	swderiv = vij*dswdr/rgrp
1050	<	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)
1041	>	if (loop .eq. PAIR_LOOP) then
1042	>	if (in_switching_region) then
1043	>	swderiv = vij*dswdr/rgrp
1044	>	fij(1) = fij(1) + swderiv*d_grp(1)
1045	>	fij(2) = fij(2) + swderiv*d_grp(2)
1046	>	fij(3) = fij(3) + swderiv*d_grp(3)
1047	>
1048	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1049	>	atom1=groupListRow(ia)
1050	>	mf = mfactRow(atom1)
1051		#ifdef IS_MPI
1052	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1053	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1054	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1052	>	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1053	>	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1054	>	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1055		#else
1056	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1057	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1058	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1056	>	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1057	>	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1058	>	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1059		#endif
1060	<	enddo
1061	<
1062	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1063	<	atom2=groupListCol(jb)
1064	<	mf = mfactCol(atom2)
1060	>	enddo
1061	>
1062	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1063	>	atom2=groupListCol(jb)
1064	>	mf = mfactCol(atom2)
1065		#ifdef IS_MPI
1066	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1067	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1068	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1066	>	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1067	>	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1068	>	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1069		#else
1070	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1071	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1072	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1070	>	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1071	>	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1072	>	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1073		#endif
1074	<	enddo
1075	<	endif
1074	>	enddo
1075	>	endif
1076
1077	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1077	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1078	>	endif
1079		endif
1080	<	end if
1080	>	endif
1081		enddo
1082	+
1083		enddo outer
1084
1085		if (update_nlist) then
#	Line 950 | Line 1139 | contains
1139
1140		if (do_pot) then
1141		! scatter/gather pot_row into the members of my column
1142	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1143	<
1142	>	do i = 1,LR_POT_TYPES
1143	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1144	>	end do
1145		! scatter/gather pot_local into all other procs
1146		! add resultant to get total pot
1147		do i = 1, nlocal
1148	<	pot_local = pot_local + pot_Temp(i)
1148	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1149	>	+ pot_Temp(1:LR_POT_TYPES,i)
1150		enddo
1151
1152		pot_Temp = 0.0_DP
1153	<
1154	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1153	>	do i = 1,LR_POT_TYPES
1154	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1155	>	end do
1156		do i = 1, nlocal
1157	<	pot_local = pot_local + pot_Temp(i)
1157	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1158	>	+ pot_Temp(1:LR_POT_TYPES,i)
1159		enddo
1160
1161		endif
1162		#endif
1163
1164	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1164	>	if (SIM_requires_postpair_calc) then
1165	>	do i = 1, nlocal
1166	>
1167	>	! we loop only over the local atoms, so we don't need row and column
1168	>	! lookups for the types
1169	>
1170	>	me_i = atid(i)
1171	>
1172	>	! is the atom electrostatic?  See if it would have an
1173	>	! electrostatic interaction with itself
1174	>	iHash = InteractionHash(me_i,me_i)
1175
1176	<	if (electrostaticSummationMethod == 3) then
974	<
1176	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1177		#ifdef IS_MPI
1178	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1179	<	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)
1178	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1179	>	t, do_pot)
1180		#else
1181	<	me_i = atid(i)
1181	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1182	>	t, do_pot)
1183		#endif
1184	<	iHash = InteractionHash(me_i,me_j)
1184	>	endif
1185	>
1186	>
1187	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1188
1189	<	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1190	<
1191	<	mu_i = getDipoleMoment(me_i)
1192	<
1193	<	!! The reaction field needs to include a self contribution
1194	<	!! to the field:
1195	<	call accumulate_self_rf(i, mu_i, eFrame)
1196	<	!! Get the reaction field contribution to the
1197	<	!! potential and torques:
1198	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1189	>	! loop over the excludes to accumulate RF stuff we've
1190	>	! left out of the normal pair loop
1191	>
1192	>	do i1 = 1, nSkipsForAtom(i)
1193	>	j = skipsForAtom(i, i1)
1194	>
1195	>	! prevent overcounting of the skips
1196	>	if (i.lt.j) then
1197	>	call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1198	>	rVal = dsqrt(ratmsq)
1199	>	call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1200		#ifdef IS_MPI
1201	<	pot_local = pot_local + rfpot
1201	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1202	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1203		#else
1204	<	pot = pot + rfpot
1205	<
1204	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1205	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1206		#endif
1207	<	endif
1208	<	enddo
1209	<	endif
1207	>	endif
1208	>	enddo
1209	>	endif
1210	>	enddo
1211		endif
1212	<
1014	<
1212	>
1213		#ifdef IS_MPI
1214	<
1214	>
1215		if (do_pot) then
1216	<	pot = pot + pot_local
1217	<	!! we assume the c code will do the allreduce to get the total potential
1020	<	!! we could do it right here if we needed to...
1216	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1217	>	mpi_comm_world,mpi_err)
1218		endif
1219	<
1219	>
1220		if (do_stress) then
1221		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1222		mpi_comm_world,mpi_err)
1223		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1224		mpi_comm_world,mpi_err)
1225		endif
1226	<
1226	>
1227		#else
1228	<
1228	>
1229		if (do_stress) then
1230		tau = tau_Temp
1231		virial = virial_Temp
1232		endif
1233	<
1233	>
1234		#endif
1235	<
1235	>
1236		end subroutine do_force_loop
1237
1238		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1239	<	eFrame, A, f, t, pot, vpair, fpair)
1239	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1240
1241	<	real( kind = dp ) :: pot, vpair, sw
1241	>	real( kind = dp ) :: vpair, sw
1242	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1243		real( kind = dp ), dimension(3) :: fpair
1244		real( kind = dp ), dimension(nLocal)   :: mfact
1245		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1052 | Line 1250 | contains
1250		logical, intent(inout) :: do_pot
1251		integer, intent(in) :: i, j
1252		real ( kind = dp ), intent(inout) :: rijsq
1253	<	real ( kind = dp )                :: r
1253	>	real ( kind = dp ), intent(inout) :: r_grp
1254		real ( kind = dp ), intent(inout) :: d(3)
1255	+	real ( kind = dp ), intent(inout) :: d_grp(3)
1256	+	real ( kind = dp ), intent(inout) :: rCut
1257	+	real ( kind = dp ) :: r
1258	+	real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1259		integer :: me_i, me_j
1260	+	integer :: k
1261
1262		integer :: iHash
1263
1264	<	r = sqrt(rijsq)
1264	>	if (useSpline) then
1265	>	call lookupUniformSpline(splineSqrt, rijsq, r)
1266	>	else
1267	>	r = sqrt(rijsq)
1268	>	endif
1269	>
1270		vpair = 0.0d0
1271		fpair(1:3) = 0.0d0
1272
#	Line 1071 | Line 1279 | contains
1279		#endif
1280
1281		iHash = InteractionHash(me_i, me_j)
1282	<
1282	>
1283		if ( iand(iHash, LJ_PAIR).ne.0 ) then
1284	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1284	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1285	>	pot(VDW_POT), f, do_pot)
1286		endif
1287	<
1287	>
1288		if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1289	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1290	<	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	<
1289	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1290	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1291		endif
1292	<
1292	>
1293		if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1294		call do_sticky_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, STICKYPOWER_PAIR).ne.0 ) then
1299		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1300	<	pot, A, f, t, do_pot)
1300	>	pot(HB_POT), A, f, t, do_pot)
1301		endif
1302	<
1302	>
1303		if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1304		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1305	<	pot, A, f, t, do_pot)
1305	>	pot(VDW_POT), A, f, t, do_pot)
1306		endif
1307
1308		if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1309	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1310	<	!           pot, A, f, t, do_pot)
1309	>	call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1310	>	pot(VDW_POT), A, f, t, do_pot)
1311		endif
1312	<
1312	>
1313		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1314	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1315	<	do_pot)
1314	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1315	>	pot(METALLIC_POT), f, do_pot)
1316		endif
1317	<
1317	>
1318		if ( iand(iHash, SHAPE_PAIR).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, SHAPE_LJ).ne.0 ) then
1324		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1325	<	pot, A, f, t, do_pot)
1325	>	pot(VDW_POT), A, f, t, do_pot)
1326		endif
1327	<
1327	>
1328	>	if ( iand(iHash, SC_PAIR).ne.0 ) then
1329	>	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1330	>	pot(METALLIC_POT), f, do_pot)
1331	>	endif
1332	>
1333		end subroutine do_pair
1334
1335	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1335	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1336		do_pot, do_stress, eFrame, A, f, t, pot)
1337
1338	<	real( kind = dp ) :: pot, sw
1338	>	real( kind = dp ) :: sw
1339	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1340		real( kind = dp ), dimension(9,nLocal) :: eFrame
1341		real (kind=dp), dimension(9,nLocal) :: A
1342		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1137 | Line 1344 | contains
1344
1345		logical, intent(inout) :: do_pot, do_stress
1346		integer, intent(in) :: i, j
1347	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1347	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1348		real ( kind = dp )                :: r, rc
1349		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1350
1351		integer :: me_i, me_j, iHash
1352
1353	<	r = sqrt(rijsq)
1353	>	if (useSpline) then
1354	>	call lookupUniformSpline(splineSqrt, rijsq, r)
1355	>	else
1356	>	r = sqrt(rijsq)
1357	>	endif
1358
1359		#ifdef IS_MPI
1360		me_i = atid_row(i)
#	Line 1156 | Line 1367 | contains
1367		iHash = InteractionHash(me_i, me_j)
1368
1369		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1370	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1370	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1371		endif
1372	+
1373	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1374	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1375	+	endif
1376
1377		end subroutine do_prepair
1378
1379
1380		subroutine do_preforce(nlocal,pot)
1381		integer :: nlocal
1382	<	real( kind = dp ) :: pot
1382	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1383
1384		if (FF_uses_EAM .and. SIM_uses_EAM) then
1385	<	call calc_EAM_preforce_Frho(nlocal,pot)
1385	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1386		endif
1387	<
1388	<
1387	>	if (FF_uses_SC .and. SIM_uses_SC) then
1388	>	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1389	>	endif
1390		end subroutine do_preforce
1391
1392
#	Line 1182 | Line 1398 | contains
1398		real( kind = dp ) :: d(3), scaled(3)
1399		integer i
1400
1401	<	d(1:3) = q_j(1:3) - q_i(1:3)
1401	>	d(1) = q_j(1) - q_i(1)
1402	>	d(2) = q_j(2) - q_i(2)
1403	>	d(3) = q_j(3) - q_i(3)
1404
1405		! Wrap back into periodic box if necessary
1406		if ( SIM_uses_PBC ) then
1407
1408		if( .not.boxIsOrthorhombic ) then
1409		! calc the scaled coordinates.
1410	+	! scaled = matmul(HmatInv, d)
1411
1412	<	scaled = matmul(HmatInv, d)
1413	<
1412	>	scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1413	>	scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1414	>	scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1415	>
1416		! wrap the scaled coordinates
1417
1418	<	scaled = scaled  - anint(scaled)
1418	>	scaled(1) = scaled(1) - dnint(scaled(1))
1419	>	scaled(2) = scaled(2) - dnint(scaled(2))
1420	>	scaled(3) = scaled(3) - dnint(scaled(3))
1421
1199	–
1422		! calc the wrapped real coordinates from the wrapped scaled
1423		! coordinates
1424	+	! d = matmul(Hmat,scaled)
1425	+	d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1426	+	d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1427	+	d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1428
1203	–	d = matmul(Hmat,scaled)
1204	–
1429		else
1430		! calc the scaled coordinates.
1431
1432	<	do i = 1, 3
1433	<	scaled(i) = d(i) * HmatInv(i,i)
1432	>	scaled(1) = d(1) * HmatInv(1,1)
1433	>	scaled(2) = d(2) * HmatInv(2,2)
1434	>	scaled(3) = d(3) * HmatInv(3,3)
1435	>
1436	>	! wrap the scaled coordinates
1437	>
1438	>	scaled(1) = scaled(1) - dnint(scaled(1))
1439	>	scaled(2) = scaled(2) - dnint(scaled(2))
1440	>	scaled(3) = scaled(3) - dnint(scaled(3))
1441
1442	<	! wrap the scaled coordinates
1442	>	! calc the wrapped real coordinates from the wrapped scaled
1443	>	! coordinates
1444
1445	<	scaled(i) = scaled(i) - anint(scaled(i))
1445	>	d(1) = scaled(1)*Hmat(1,1)
1446	>	d(2) = scaled(2)*Hmat(2,2)
1447	>	d(3) = scaled(3)*Hmat(3,3)
1448
1215	–	! calc the wrapped real coordinates from the wrapped scaled
1216	–	! coordinates
1217	–
1218	–	d(i) = scaled(i)*Hmat(i,i)
1219	–	enddo
1449		endif
1450
1451		endif
1452
1453	<	r_sq = dot_product(d,d)
1453	>	r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1454
1455		end subroutine get_interatomic_vector
1456
#	Line 1253 | Line 1482 | contains
1482		pot_Col = 0.0_dp
1483		pot_Temp = 0.0_dp
1484
1256	–	rf_Row = 0.0_dp
1257	–	rf_Col = 0.0_dp
1258	–	rf_Temp = 0.0_dp
1259	–
1485		#endif
1486
1487		if (FF_uses_EAM .and. SIM_uses_EAM) then
1488		call clean_EAM()
1489		endif
1490
1266	–	rf = 0.0_dp
1491		tau_Temp = 0.0_dp
1492		virial_Temp = 0.0_dp
1493		end subroutine zero_work_arrays
#	Line 1357 | Line 1581 | contains
1581
1582		function FF_RequiresPrepairCalc() result(doesit)
1583		logical :: doesit
1584	<	doesit = FF_uses_EAM
1584	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1585	>	.or. FF_uses_MEAM
1586		end function FF_RequiresPrepairCalc
1587
1363	–	function FF_RequiresPostpairCalc() result(doesit)
1364	–	logical :: doesit
1365	–	if (electrostaticSummationMethod == 3) doesit = .true.
1366	–	end function FF_RequiresPostpairCalc
1367	–
1588		#ifdef PROFILE
1589		function getforcetime() result(totalforcetime)
1590		real(kind=dp) :: totalforcetime
#	Line 1397 | Line 1617 | end module doForces
1617
1618		end subroutine add_stress_tensor
1619
1620	+	subroutine setupSplineSqrt(rmax)
1621	+	real(kind=dp), intent(in) :: rmax
1622	+	real(kind=dp), dimension(np) :: xvals, yvals
1623	+	real(kind=dp) :: r2_1, r2_n, dx, r2
1624	+	integer :: i
1625	+
1626	+	r2_1 = 0.5d0
1627	+	r2_n = rmax*rmax
1628	+
1629	+	dx = (r2_n-r2_1) / dble(np-1)
1630	+
1631	+	do i = 1, np
1632	+	r2 = r2_1 + dble(i-1)*dx
1633	+	xvals(i) = r2
1634	+	yvals(i) = dsqrt(r2)
1635	+	enddo
1636	+
1637	+	call newSpline(splineSqrt, xvals, yvals, .true.)
1638	+
1639	+	haveSplineSqrt = .true.
1640	+	return
1641	+	end subroutine setupSplineSqrt
1642	+
1643	+	subroutine deleteSplineSqrt()
1644	+	call deleteSpline(splineSqrt)
1645	+	haveSplineSqrt = .false.
1646	+	return
1647	+	end subroutine deleteSplineSqrt
1648	+
1649		end module doForces

Diff Legend

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