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 2284 by gezelter, Wed Sep 7 19:18:54 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.37 2005-09-07 19:18:54 gezelter Exp $, $Date: 2005-09-07 19:18:54 $, $Name: not supported by cvs2svn $, $Revision: 1.37 $
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
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_RF
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_RF
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 :: corrMethod
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
108	–	public :: createInteractionHash
109	–	public :: createGtypeCutoffMap
110	–	public :: getStickyCut
111	–	public :: getStickyPowerCut
112	–	public :: getGayBerneCut
113	–	public :: getEAMCut
114	–	public :: getShapeCut
126
127		#ifdef PROFILE
128		public :: getforcetime
#	Line 124 | 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 134 | 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
144	–	integer, intent(out) :: status
162		integer :: i
163		integer :: j
164		integer :: iHash
#	Line 153 | 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 160 | 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
165	–	status = 0
166	–
186		if (.not. associated(atypes)) then
187	<	call handleError("atype", "atypes was not present before call of createInteractionHash!")
169	<	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 193 | 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 206 | 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 227 | 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 246 | Line 278 | contains
278		haveInteractionHash = .true.
279		end subroutine createInteractionHash
280
281	<	subroutine createGtypeCutoffMap(stat)
281	>	subroutine createGtypeCutoffMap()
282
251	–	integer, intent(out), optional :: stat
283		logical :: i_is_LJ
284		logical :: i_is_Elect
285		logical :: i_is_Sticky
#	Line 256 | 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
295	<	real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
294	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
295	>	integer :: nGroupsInRow
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	<
304	>	#ifdef IS_MPI
305	>	nGroupsInRow = getNgroupsInRow(plan_group_row)
306	>	nGroupsInCol = getNgroupsInCol(plan_group_col)
307	>	#endif
308		nAtypes = getSize(atypes)
309	<
309	>	! Set all of the initial cutoffs to zero.
310	>	atypeMaxCutoff = 0.0_dp
311		do i = 1, nAtypes
312		if (SimHasAtype(i)) then
313		call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
#	Line 283 | 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	<
321	<	atypeMaxCutoff(i) = 0.0_dp
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
320	>	call getElementProperty(atypes, i, "is_SC", i_is_SC)
321	>
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
322	–
323	–	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))
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 351 | 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)
448	<	endif
446	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
447	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
448	>	endif
449		enddo
450	<	if (nGroupTypes.eq.0) then
451	<	nGroupTypes = nGroupTypes + 1
452	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
453	<	groupToGtype(i) = nGroupTypes
450	>	if (nGroupTypesRow.eq.0) then
451	>	nGroupTypesRow = nGroupTypesRow + 1
452	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
453	>	groupToGtypeRow(i) = nGroupTypesRow
454		else
455	<	do g = 1, nGroupTypes
456	<	if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
457	<	nGroupTypes = nGroupTypes + 1
458	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
459	<	groupToGtype(i) = nGroupTypes
368	<	else
369	<	groupToGtype(i) = g
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	<
468	>	enddo
469	>
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, 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	>	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
396	–	skin = defaultRlist - defaultRcut
397	–	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	<
404	<	end subroutine createGtypeCutoffMap
405	<
406	<	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
407	<	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
408	<	integer, intent(in) :: cutPolicy
409	<
410	<	defaultRcut = defRcut
411	<	defaultRsw = defRsw
412	<	defaultRlist = defRlist
413	<	cutoffPolicy = cutPolicy
414	<	end subroutine setDefaultCutoffs
415	<
416	<	subroutine setCutoffPolicy(cutPolicy)
565	>	end subroutine createGtypeCutoffMap
566
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	+
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	+	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	<
616	>	haveCutoffPolicy = .true.
617	>	haveGtypeCutoffMap = .false.
618	>
619		end subroutine setCutoffPolicy
620	+
621	+	subroutine setElectrostaticMethod( thisESM )
622	+
623	+	integer, intent(in) :: thisESM
624	+
625	+	electrostaticSummationMethod = thisESM
626	+	haveElectrostaticSummationMethod = .true.
627
628	+	end subroutine setElectrostaticMethod
629	+
630	+	subroutine setSkinThickness( thisSkin )
631
632	<	subroutine setSimVariables()
633	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
634	<	SIM_uses_EAM = SimUsesEAM()
635	<	SIM_uses_RF = SimUsesRF()
636	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
637	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
638	<	SIM_uses_PBC = SimUsesPBC()
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
433	–	haveSIMvariables = .true.
434	–
435	–	return
436	–	end subroutine setSimVariables
437	–
653		subroutine doReadyCheck(error)
654		integer, intent(out) :: error
440	–
655		integer :: myStatus
656
657		error = 0
658
659		if (.not. haveInteractionHash) then
660	<	myStatus = 0
447	<	call createInteractionHash(myStatus)
448	<	if (myStatus .ne. 0) then
449	<	write(default_error, *) 'createInteractionHash failed in doForces!'
450	<	error = -1
451	<	return
452	<	endif
660	>	call createInteractionHash()
661		endif
662
663		if (.not. haveGtypeCutoffMap) then
664	<	myStatus = 0
457	<	call createGtypeCutoffMap(myStatus)
458	<	if (myStatus .ne. 0) then
459	<	write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
460	<	error = -1
461	<	return
462	<	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
469	–	!  if (.not. haveRlist) then
470	–	!     write(default_error, *) 'rList has not been set in doForces!'
471	–	!     error = -1
472	–	!     return
473	–	!  endif
474	–
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 495 | Line 704 | contains
704		end subroutine doReadyCheck
705
706
707	<	subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
707	>	subroutine init_FF(thisStat)
708
500	–	logical, intent(in) :: use_RF
501	–	logical, intent(in) :: use_UW
502	–	logical, intent(in) :: use_DW
709		integer, intent(out) :: thisStat
710		integer :: my_status, nMatches
505	–	integer :: corrMethod
711		integer, pointer :: MatchList(:) => null()
507	–	real(kind=dp) :: rcut, rrf, rt, dielect
712
713		!! assume things are copacetic, unless they aren't
714		thisStat = 0
715
512	–	!! Fortran's version of a cast:
513	–	FF_uses_RF = use_RF
514	–
515	–	!! set the electrostatic correction method
516	–	if (use_UW) then
517	–	corrMethod = 1
518	–	elseif (use_DW) then
519	–	corrMethod = 2
520	–	else
521	–	corrMethod = 0
522	–	endif
523	–
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 531 | 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 547 | 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
551	–	haveSaneForceField = .true.
746
747	<	!! check to make sure the FF_uses_RF setting makes sense
554	<
555	<	if (FF_uses_Dipoles) then
556	<	if (FF_uses_RF) then
557	<	dielect = getDielect()
558	<	call initialize_rf(dielect)
559	<	endif
560	<	else
561	<	if (FF_uses_RF) then
562	<	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
563	<	thisStat = -1
564	<	haveSaneForceField = .false.
565	<	return
566	<	endif
567	<	endif
747	>	haveSaneForceField = .true.
748
749		if (FF_uses_EAM) then
750		call init_EAM_FF(my_status)
#	Line 576 | Line 756 | contains
756		end if
757		endif
758
579	–	if (FF_uses_GayBerne) then
580	–	call check_gb_pair_FF(my_status)
581	–	if (my_status .ne. 0) then
582	–	thisStat = -1
583	–	haveSaneForceField = .false.
584	–	return
585	–	endif
586	–	endif
587	–
759		if (.not. haveNeighborList) then
760		!! Create neighbor lists
761		call expandNeighborList(nLocal, my_status)
#	Line 618 | 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 639 | 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 649 | 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 713 | 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 774 | 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 796 | 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)
864	<	fij(2) = fij(2) + swderiv*d_grp(2)
865	<	fij(3) = fij(3) + swderiv*d_grp(3)
866	<
867	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
868	<	atom1=groupListRow(ia)
869	<	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 956 | 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 (FF_uses_RF .and. SIM_uses_RF) then
980	<
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)
984	<	do i = 1,nlocal
985	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
986	<	end do
987	<	#endif
988	<
989	<	do i = 1, nLocal
990	<
991	<	rfpot = 0.0_DP
992	<	#ifdef IS_MPI
993	<	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	<
1020	<
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
1026	<	!! 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 1058 | 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 ) :: ebalance
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 1078 | 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, corrMethod)
1089	<
1090	<	if (FF_uses_RF .and. SIM_uses_RF) then
1091	<
1092	<	! CHECK ME (RF needs to know about all electrostatic types)
1093	<	call accumulate_rf(i, j, r, eFrame, sw)
1094	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1095	<	endif
1096	<
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 1144 | 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	+	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)
1361		me_j = atid_col(j)
#	Line 1161 | 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
#	Line 1169 | Line 1379 | contains
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 1187 | 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
1204	–
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
1208	–	d = matmul(Hmat,scaled)
1209	–
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))
1446	<
1447	<	! calc the wrapped real coordinates from the wrapped scaled
1221	<	! coordinates
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
1223	–	d(i) = scaled(i)*Hmat(i,i)
1224	–	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 1258 | Line 1482 | contains
1482		pot_Col = 0.0_dp
1483		pot_Temp = 0.0_dp
1484
1261	–	rf_Row = 0.0_dp
1262	–	rf_Col = 0.0_dp
1263	–	rf_Temp = 0.0_dp
1264	–
1485		#endif
1486
1487		if (FF_uses_EAM .and. SIM_uses_EAM) then
1488		call clean_EAM()
1489		endif
1490
1271	–	rf = 0.0_dp
1491		tau_Temp = 0.0_dp
1492		virial_Temp = 0.0_dp
1493		end subroutine zero_work_arrays
#	Line 1362 | 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
1368	–	function FF_RequiresPostpairCalc() result(doesit)
1369	–	logical :: doesit
1370	–	doesit = FF_uses_RF
1371	–	end function FF_RequiresPostpairCalc
1372	–
1588		#ifdef PROFILE
1589		function getforcetime() result(totalforcetime)
1590		real(kind=dp) :: totalforcetime
#	Line 1402 | 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