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 2285 by gezelter, Wed Sep 7 20:46:46 2005 UTC vs.
Revision 2917 by chrisfen, Mon Jul 3 13:18:43 2006 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.38 2005-09-07 20:46:39 gezelter Exp $, $Date: 2005-09-07 20:46:39 $, $Name: not supported by cvs2svn $, $Revision: 1.38 $
48	>	!! @version $Id: doForces.F90,v 1.84 2006-07-03 13:18:43 chrisfen Exp $, $Date: 2006-07-03 13:18:43 $, $Name: not supported by cvs2svn $, $Revision: 1.84 $
49
50
51		module doForces
#	Line 58 | Line 58 | module doForces
58		use lj
59		use sticky
60		use electrostatic_module
61	<	use reaction_field
62	<	use gb_pair
61	>	use gayberne
62		use shapes
63		use vector_class
64		use eam
65	+	use suttonchen
66		use status
67		#ifdef IS_MPI
68		use mpiSimulation
#	Line 72 | Line 72 | module doForces
72		PRIVATE
73
74		#define __FORTRAN90
75	–	#include "UseTheForce/fSwitchingFunction.h"
75		#include "UseTheForce/fCutoffPolicy.h"
76		#include "UseTheForce/DarkSide/fInteractionMap.h"
77	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78
79	–
79		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
80		INTEGER, PARAMETER:: PAIR_LOOP    = 2
81
#	Line 85 | Line 84 | module doForces
84		logical, save :: haveSaneForceField = .false.
85		logical, save :: haveInteractionHash = .false.
86		logical, save :: haveGtypeCutoffMap = .false.
87	<	logical, save :: haveRlist = .false.
87	>	logical, save :: haveDefaultCutoffs = .false.
88	>	logical, save :: haveSkinThickness = .false.
89	>	logical, save :: haveElectrostaticSummationMethod = .false.
90	>	logical, save :: haveCutoffPolicy = .false.
91	>	logical, save :: VisitCutoffsAfterComputing = .false.
92	>	logical, save :: do_box_dipole = .false.
93
94		logical, save :: FF_uses_DirectionalAtoms
95		logical, save :: FF_uses_Dipoles
96		logical, save :: FF_uses_GayBerne
97		logical, save :: FF_uses_EAM
98	<	logical, save :: FF_uses_RF
98	>	logical, save :: FF_uses_SC
99	>	logical, save :: FF_uses_MEAM
100	>
101
102		logical, save :: SIM_uses_DirectionalAtoms
103		logical, save :: SIM_uses_EAM
104	<	logical, save :: SIM_uses_RF
104	>	logical, save :: SIM_uses_SC
105	>	logical, save :: SIM_uses_MEAM
106		logical, save :: SIM_requires_postpair_calc
107		logical, save :: SIM_requires_prepair_calc
108		logical, save :: SIM_uses_PBC
109
110	<	integer, save :: corrMethod
110	>	integer, save :: electrostaticSummationMethod
111	>	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
112
113	+	real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
114	+	real(kind=dp), save :: skinThickness
115	+	logical, save :: defaultDoShift
116	+
117		public :: init_FF
118	<	public :: setDefaultCutoffs
118	>	public :: setCutoffs
119	>	public :: cWasLame
120	>	public :: setElectrostaticMethod
121	>	public :: setBoxDipole
122	>	public :: getBoxDipole
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
139	<
153	>	real(kind=dp), dimension(3) :: boxDipole
154	>
155		contains
156
157	<	subroutine createInteractionHash(status)
157	>	subroutine createInteractionHash()
158		integer :: nAtypes
144	–	integer, intent(out) :: status
159		integer :: i
160		integer :: j
161		integer :: iHash
#	Line 153 | Line 167 | contains
167		logical :: i_is_GB
168		logical :: i_is_EAM
169		logical :: i_is_Shape
170	+	logical :: i_is_SC
171	+	logical :: i_is_MEAM
172		logical :: j_is_LJ
173		logical :: j_is_Elect
174		logical :: j_is_Sticky
#	Line 160 | Line 176 | contains
176		logical :: j_is_GB
177		logical :: j_is_EAM
178		logical :: j_is_Shape
179	+	logical :: j_is_SC
180	+	logical :: j_is_MEAM
181		real(kind=dp) :: myRcut
182
165	–	status = 0
166	–
183		if (.not. associated(atypes)) then
184	<	call handleError("atype", "atypes was not present before call of createInteractionHash!")
169	<	status = -1
184	>	call handleError("doForces", "atypes was not present before call of createInteractionHash!")
185		return
186		endif
187
188		nAtypes = getSize(atypes)
189
190		if (nAtypes == 0) then
191	<	status = -1
191	>	call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
192		return
193		end if
194
195		if (.not. allocated(InteractionHash)) then
196		allocate(InteractionHash(nAtypes,nAtypes))
197	+	else
198	+	deallocate(InteractionHash)
199	+	allocate(InteractionHash(nAtypes,nAtypes))
200		endif
201
202		if (.not. allocated(atypeMaxCutoff)) then
203		allocate(atypeMaxCutoff(nAtypes))
204	+	else
205	+	deallocate(atypeMaxCutoff)
206	+	allocate(atypeMaxCutoff(nAtypes))
207		endif
208
209		do i = 1, nAtypes
#	Line 193 | Line 214 | contains
214		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
215		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
216		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
217	+	call getElementProperty(atypes, i, "is_SC", i_is_SC)
218	+	call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
219
220		do j = i, nAtypes
221
#	Line 206 | Line 229 | contains
229		call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
230		call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
231		call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
232	+	call getElementProperty(atypes, j, "is_SC", j_is_SC)
233	+	call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
234
235		if (i_is_LJ .and. j_is_LJ) then
236		iHash = ior(iHash, LJ_PAIR)
#	Line 225 | Line 250 | contains
250
251		if (i_is_EAM .and. j_is_EAM) then
252		iHash = ior(iHash, EAM_PAIR)
253	+	endif
254	+
255	+	if (i_is_SC .and. j_is_SC) then
256	+	iHash = ior(iHash, SC_PAIR)
257		endif
258
259		if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
#	Line 246 | Line 275 | contains
275		haveInteractionHash = .true.
276		end subroutine createInteractionHash
277
278	<	subroutine createGtypeCutoffMap(stat)
278	>	subroutine createGtypeCutoffMap()
279
251	–	integer, intent(out), optional :: stat
280		logical :: i_is_LJ
281		logical :: i_is_Elect
282		logical :: i_is_Sticky
#	Line 256 | Line 284 | contains
284		logical :: i_is_GB
285		logical :: i_is_EAM
286		logical :: i_is_Shape
287	+	logical :: i_is_SC
288	+	logical :: GtypeFound
289
290		integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
291	<	integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
292	<	real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
291	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
292	>	integer :: nGroupsInRow
293	>	integer :: nGroupsInCol
294	>	integer :: nGroupTypesRow,nGroupTypesCol
295	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
296		real(kind=dp) :: biggestAtypeCutoff
297
265	–	stat = 0
298		if (.not. haveInteractionHash) then
299	<	call createInteractionHash(myStatus)
268	<	if (myStatus .ne. 0) then
269	<	write(default_error, *) 'createInteractionHash failed in doForces!'
270	<	stat = -1
271	<	return
272	<	endif
299	>	call createInteractionHash()
300		endif
301	<
301	>	#ifdef IS_MPI
302	>	nGroupsInRow = getNgroupsInRow(plan_group_row)
303	>	nGroupsInCol = getNgroupsInCol(plan_group_col)
304	>	#endif
305		nAtypes = getSize(atypes)
306	<
306	>	! Set all of the initial cutoffs to zero.
307	>	atypeMaxCutoff = 0.0_dp
308		do i = 1, nAtypes
309		if (SimHasAtype(i)) then
310		call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
#	Line 283 | Line 314 | contains
314		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
315		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
316		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
317	<
318	<	atypeMaxCutoff(i) = 0.0_dp
319	<	if (i_is_LJ) then
320	<	thisRcut = getSigma(i) * 2.5_dp
321	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322	<	endif
323	<	if (i_is_Elect) then
324	<	thisRcut = defaultRcut
325	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326	<	endif
327	<	if (i_is_Sticky) then
328	<	thisRcut = getStickyCut(i)
329	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330	<	endif
331	<	if (i_is_StickyP) then
332	<	thisRcut = getStickyPowerCut(i)
333	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334	<	endif
335	<	if (i_is_GB) then
336	<	thisRcut = getGayBerneCut(i)
337	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
338	<	endif
339	<	if (i_is_EAM) then
340	<	thisRcut = getEAMCut(i)
341	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
342	<	endif
343	<	if (i_is_Shape) then
344	<	thisRcut = getShapeCut(i)
345	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317	>	call getElementProperty(atypes, i, "is_SC", i_is_SC)
318	>
319	>	if (haveDefaultCutoffs) then
320	>	atypeMaxCutoff(i) = defaultRcut
321	>	else
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
349	>	endif
350	>	if (i_is_SC) then
351	>	thisRcut = getSCCut(i)
352	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
353	>	endif
354		endif
355	<
355	>
356		if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
357		biggestAtypeCutoff = atypeMaxCutoff(i)
358		endif
359	+
360		endif
361		enddo
322	–
323	–	nGroupTypes = 0
362
363		istart = 1
364	+	jstart = 1
365		#ifdef IS_MPI
366		iend = nGroupsInRow
367	+	jend = nGroupsInCol
368		#else
369		iend = nGroups
370	+	jend = nGroups
371		#endif
372
373		!! allocate the groupToGtype and gtypeMaxCutoff here.
374	<	if(.not.allocated(groupToGtype)) then
375	<	allocate(groupToGtype(iend))
376	<	allocate(groupMaxCutoff(iend))
377	<	allocate(gtypeMaxCutoff(iend))
374	>	if(.not.allocated(groupToGtypeRow)) then
375	>	!  allocate(groupToGtype(iend))
376	>	allocate(groupToGtypeRow(iend))
377	>	else
378	>	deallocate(groupToGtypeRow)
379	>	allocate(groupToGtypeRow(iend))
380		endif
381	+	if(.not.allocated(groupMaxCutoffRow)) then
382	+	allocate(groupMaxCutoffRow(iend))
383	+	else
384	+	deallocate(groupMaxCutoffRow)
385	+	allocate(groupMaxCutoffRow(iend))
386	+	end if
387	+
388	+	if(.not.allocated(gtypeMaxCutoffRow)) then
389	+	allocate(gtypeMaxCutoffRow(iend))
390	+	else
391	+	deallocate(gtypeMaxCutoffRow)
392	+	allocate(gtypeMaxCutoffRow(iend))
393	+	endif
394	+
395	+
396	+	#ifdef IS_MPI
397	+	! We only allocate new storage if we are in MPI because Ncol /= Nrow
398	+	if(.not.associated(groupToGtypeCol)) then
399	+	allocate(groupToGtypeCol(jend))
400	+	else
401	+	deallocate(groupToGtypeCol)
402	+	allocate(groupToGtypeCol(jend))
403	+	end if
404	+
405	+	if(.not.associated(groupMaxCutoffCol)) then
406	+	allocate(groupMaxCutoffCol(jend))
407	+	else
408	+	deallocate(groupMaxCutoffCol)
409	+	allocate(groupMaxCutoffCol(jend))
410	+	end if
411	+	if(.not.associated(gtypeMaxCutoffCol)) then
412	+	allocate(gtypeMaxCutoffCol(jend))
413	+	else
414	+	deallocate(gtypeMaxCutoffCol)
415	+	allocate(gtypeMaxCutoffCol(jend))
416	+	end if
417	+
418	+	groupMaxCutoffCol = 0.0_dp
419	+	gtypeMaxCutoffCol = 0.0_dp
420	+
421	+	#endif
422	+	groupMaxCutoffRow = 0.0_dp
423	+	gtypeMaxCutoffRow = 0.0_dp
424	+
425	+
426		!! first we do a single loop over the cutoff groups to find the
427		!! largest cutoff for any atypes present in this group.  We also
428		!! create gtypes at this point.
429
430	<	tol = 1.0d-6
431	<
430	>	tol = 1.0e-6_dp
431	>	nGroupTypesRow = 0
432	>	nGroupTypesCol = 0
433		do i = istart, iend
434		n_in_i = groupStartRow(i+1) - groupStartRow(i)
435	<	groupMaxCutoff(i) = 0.0_dp
435	>	groupMaxCutoffRow(i) = 0.0_dp
436		do ia = groupStartRow(i), groupStartRow(i+1)-1
437		atom1 = groupListRow(ia)
438		#ifdef IS_MPI
#	Line 351 | Line 440 | contains
440		#else
441		me_i = atid(atom1)
442		#endif
443	<	if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
444	<	groupMaxCutoff(i)=atypeMaxCutoff(me_i)
443	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
444	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
445	>	endif
446	>	enddo
447	>	if (nGroupTypesRow.eq.0) then
448	>	nGroupTypesRow = nGroupTypesRow + 1
449	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
450	>	groupToGtypeRow(i) = nGroupTypesRow
451	>	else
452	>	GtypeFound = .false.
453	>	do g = 1, nGroupTypesRow
454	>	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
455	>	groupToGtypeRow(i) = g
456	>	GtypeFound = .true.
457	>	endif
458	>	enddo
459	>	if (.not.GtypeFound) then
460	>	nGroupTypesRow = nGroupTypesRow + 1
461	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
462	>	groupToGtypeRow(i) = nGroupTypesRow
463		endif
464	+	endif
465	+	enddo
466	+
467	+	#ifdef IS_MPI
468	+	do j = jstart, jend
469	+	n_in_j = groupStartCol(j+1) - groupStartCol(j)
470	+	groupMaxCutoffCol(j) = 0.0_dp
471	+	do ja = groupStartCol(j), groupStartCol(j+1)-1
472	+	atom1 = groupListCol(ja)
473	+
474	+	me_j = atid_col(atom1)
475	+
476	+	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
477	+	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
478	+	endif
479		enddo
480	<	if (nGroupTypes.eq.0) then
481	<	nGroupTypes = nGroupTypes + 1
482	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
483	<	groupToGtype(i) = nGroupTypes
480	>
481	>	if (nGroupTypesCol.eq.0) then
482	>	nGroupTypesCol = nGroupTypesCol + 1
483	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
484	>	groupToGtypeCol(j) = nGroupTypesCol
485		else
486	<	do g = 1, nGroupTypes
487	<	if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
488	<	nGroupTypes = nGroupTypes + 1
489	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
490	<	groupToGtype(i) = nGroupTypes
368	<	else
369	<	groupToGtype(i) = g
486	>	GtypeFound = .false.
487	>	do g = 1, nGroupTypesCol
488	>	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
489	>	groupToGtypeCol(j) = g
490	>	GtypeFound = .true.
491		endif
492		enddo
493	+	if (.not.GtypeFound) then
494	+	nGroupTypesCol = nGroupTypesCol + 1
495	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
496	+	groupToGtypeCol(j) = nGroupTypesCol
497	+	endif
498		endif
499	<	enddo
500	<
499	>	enddo
500	>
501	>	#else
502	>	! Set pointers to information we just found
503	>	nGroupTypesCol = nGroupTypesRow
504	>	groupToGtypeCol => groupToGtypeRow
505	>	gtypeMaxCutoffCol => gtypeMaxCutoffRow
506	>	groupMaxCutoffCol => groupMaxCutoffRow
507	>	#endif
508	>
509		!! allocate the gtypeCutoffMap here.
510	<	allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
510	>	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
511		!! then we do a double loop over all the group TYPES to find the cutoff
512		!! map between groups of two types
513	<
514	<	do i = 1, nGroupTypes
515	<	do j = 1, nGroupTypes
513	>	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
514	>
515	>	do i = 1, nGroupTypesRow
516	>	do j = 1, nGroupTypesCol
517
383	–	write(*,*) 'cutoffPolicy = ', cutoffPolicy
518		select case(cutoffPolicy)
519		case(TRADITIONAL_CUTOFF_POLICY)
520	<	thisRcut = maxval(gtypeMaxCutoff)
520	>	thisRcut = tradRcut
521		case(MIX_CUTOFF_POLICY)
522	<	thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
522	>	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
523		case(MAX_CUTOFF_POLICY)
524	<	thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
524	>	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
525		case default
526		call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
527		return
528		end select
529		gtypeCutoffMap(i,j)%rcut = thisRcut
530	+
531	+	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
532	+
533		gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
397	–	skin = defaultRlist - defaultRcut
398	–	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
534
535	+	if (.not.haveSkinThickness) then
536	+	skinThickness = 1.0_dp
537	+	endif
538	+
539	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
540	+
541	+	! sanity check
542	+
543	+	if (haveDefaultCutoffs) then
544	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
545	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
546	+	endif
547	+	endif
548		enddo
549		enddo
550	+
551	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
552	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
553	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
554	+	#ifdef IS_MPI
555	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
556	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
557	+	#endif
558	+	groupMaxCutoffCol => null()
559	+	gtypeMaxCutoffCol => null()
560
561		haveGtypeCutoffMap = .true.
562	<
405	<	end subroutine createGtypeCutoffMap
406	<
407	<	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
408	<	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
409	<	integer, intent(in) :: cutPolicy
410	<
411	<	defaultRcut = defRcut
412	<	defaultRsw = defRsw
413	<	defaultRlist = defRlist
414	<	cutoffPolicy = cutPolicy
415	<	end subroutine setDefaultCutoffs
416	<
417	<	subroutine setCutoffPolicy(cutPolicy)
562	>	end subroutine createGtypeCutoffMap
563
564	+	subroutine setCutoffs(defRcut, defRsw)
565	+
566	+	real(kind=dp),intent(in) :: defRcut, defRsw
567	+	character(len = statusMsgSize) :: errMsg
568	+	integer :: localError
569	+
570	+	defaultRcut = defRcut
571	+	defaultRsw = defRsw
572	+
573	+	defaultDoShift = .false.
574	+	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
575	+
576	+	write(errMsg, *) &
577	+	'cutoffRadius and switchingRadius are set to the same', newline &
578	+	// tab, 'value.  OOPSE will use shifted ', newline &
579	+	// tab, 'potentials instead of switching functions.'
580	+
581	+	call handleInfo("setCutoffs", errMsg)
582	+
583	+	defaultDoShift = .true.
584	+
585	+	endif
586	+
587	+	localError = 0
588	+	call setLJDefaultCutoff( defaultRcut, defaultDoShift )
589	+	call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
590	+	call setCutoffEAM( defaultRcut )
591	+	call setCutoffSC( defaultRcut )
592	+	call set_switch(defaultRsw, defaultRcut)
593	+	call setHmatDangerousRcutValue(defaultRcut)
594	+
595	+	haveDefaultCutoffs = .true.
596	+	haveGtypeCutoffMap = .false.
597	+
598	+	end subroutine setCutoffs
599	+
600	+	subroutine cWasLame()
601	+
602	+	VisitCutoffsAfterComputing = .true.
603	+	return
604	+
605	+	end subroutine cWasLame
606	+
607	+	subroutine setCutoffPolicy(cutPolicy)
608	+
609		integer, intent(in) :: cutPolicy
610	+
611		cutoffPolicy = cutPolicy
612	<	call createGtypeCutoffMap()
613	<
612	>	haveCutoffPolicy = .true.
613	>	haveGtypeCutoffMap = .false.
614	>
615		end subroutine setCutoffPolicy
616	+
617	+	subroutine setBoxDipole()
618	+
619	+	do_box_dipole = .true.
620
621	<
426	<	subroutine setSimVariables()
427	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
428	<	SIM_uses_EAM = SimUsesEAM()
429	<	SIM_uses_RF = SimUsesRF()
430	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
431	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
432	<	SIM_uses_PBC = SimUsesPBC()
621	>	end subroutine setBoxDipole
622
623	<	haveSIMvariables = .true.
623	>	subroutine getBoxDipole( box_dipole )
624
625	<	return
626	<	end subroutine setSimVariables
625	>	real(kind=dp), intent(inout), dimension(3) :: box_dipole
626	>
627	>	box_dipole = boxDipole
628	>
629	>	end subroutine getBoxDipole
630	>
631	>	subroutine setElectrostaticMethod( thisESM )
632	>
633	>	integer, intent(in) :: thisESM
634	>
635	>	electrostaticSummationMethod = thisESM
636	>	haveElectrostaticSummationMethod = .true.
637	>
638	>	end subroutine setElectrostaticMethod
639	>
640	>	subroutine setSkinThickness( thisSkin )
641	>
642	>	real(kind=dp), intent(in) :: thisSkin
643	>
644	>	skinThickness = thisSkin
645	>	haveSkinThickness = .true.
646	>	haveGtypeCutoffMap = .false.
647	>
648	>	end subroutine setSkinThickness
649	>
650	>	subroutine setSimVariables()
651	>	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
652	>	SIM_uses_EAM = SimUsesEAM()
653	>	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
654	>	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
655	>	SIM_uses_PBC = SimUsesPBC()
656	>	SIM_uses_SC = SimUsesSC()
657
658	+	haveSIMvariables = .true.
659	+
660	+	return
661	+	end subroutine setSimVariables
662	+
663		subroutine doReadyCheck(error)
664		integer, intent(out) :: error
441	–
665		integer :: myStatus
666
667		error = 0
668
669		if (.not. haveInteractionHash) then
670	<	myStatus = 0
448	<	call createInteractionHash(myStatus)
449	<	if (myStatus .ne. 0) then
450	<	write(default_error, *) 'createInteractionHash failed in doForces!'
451	<	error = -1
452	<	return
453	<	endif
670	>	call createInteractionHash()
671		endif
672
673		if (.not. haveGtypeCutoffMap) then
674	<	myStatus = 0
458	<	call createGtypeCutoffMap(myStatus)
459	<	if (myStatus .ne. 0) then
460	<	write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
461	<	error = -1
462	<	return
463	<	endif
674	>	call createGtypeCutoffMap()
675		endif
676
677	+	if (VisitCutoffsAfterComputing) then
678	+	call set_switch(largestRcut, largestRcut)
679	+	call setHmatDangerousRcutValue(largestRcut)
680	+	call setCutoffEAM(largestRcut)
681	+	call setCutoffSC(largestRcut)
682	+	VisitCutoffsAfterComputing = .false.
683	+	endif
684	+
685		if (.not. haveSIMvariables) then
686		call setSimVariables()
687		endif
688
470	–	!  if (.not. haveRlist) then
471	–	!     write(default_error, *) 'rList has not been set in doForces!'
472	–	!     error = -1
473	–	!     return
474	–	!  endif
475	–
689		if (.not. haveNeighborList) then
690		write(default_error, *) 'neighbor list has not been initialized in doForces!'
691		error = -1
692		return
693		end if
694	<
694	>
695		if (.not. haveSaneForceField) then
696		write(default_error, *) 'Force Field is not sane in doForces!'
697		error = -1
698		return
699		end if
700	<
700	>
701		#ifdef IS_MPI
702		if (.not. isMPISimSet()) then
703		write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
#	Line 496 | Line 709 | contains
709		end subroutine doReadyCheck
710
711
712	<	subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
712	>	subroutine init_FF(thisStat)
713
501	–	logical, intent(in) :: use_RF
502	–	logical, intent(in) :: use_UW
503	–	logical, intent(in) :: use_DW
714		integer, intent(out) :: thisStat
715		integer :: my_status, nMatches
506	–	integer :: corrMethod
716		integer, pointer :: MatchList(:) => null()
508	–	real(kind=dp) :: rcut, rrf, rt, dielect
717
718		!! assume things are copacetic, unless they aren't
719		thisStat = 0
720
513	–	!! Fortran's version of a cast:
514	–	FF_uses_RF = use_RF
515	–
516	–	!! set the electrostatic correction method
517	–	if (use_UW) then
518	–	corrMethod = 1
519	–	elseif (use_DW) then
520	–	corrMethod = 2
521	–	else
522	–	corrMethod = 0
523	–	endif
524	–
721		!! init_FF is called *after* all of the atom types have been
722		!! defined in atype_module using the new_atype subroutine.
723		!!
#	Line 532 | Line 728 | contains
728		FF_uses_Dipoles = .false.
729		FF_uses_GayBerne = .false.
730		FF_uses_EAM = .false.
731	+	FF_uses_SC = .false.
732
733		call getMatchingElementList(atypes, "is_Directional", .true., &
734		nMatches, MatchList)
#	Line 548 | Line 745 | contains
745		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
746		if (nMatches .gt. 0) FF_uses_EAM = .true.
747
748	+	call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
749	+	if (nMatches .gt. 0) FF_uses_SC = .true.
750
552	–	haveSaneForceField = .true.
751
752	<	!! check to make sure the FF_uses_RF setting makes sense
555	<
556	<	if (FF_uses_Dipoles) then
557	<	if (FF_uses_RF) then
558	<	dielect = getDielect()
559	<	call initialize_rf(dielect)
560	<	endif
561	<	else
562	<	if (FF_uses_RF) then
563	<	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
564	<	thisStat = -1
565	<	haveSaneForceField = .false.
566	<	return
567	<	endif
568	<	endif
752	>	haveSaneForceField = .true.
753
754		if (FF_uses_EAM) then
755		call init_EAM_FF(my_status)
#	Line 577 | Line 761 | contains
761		end if
762		endif
763
580	–	if (FF_uses_GayBerne) then
581	–	call check_gb_pair_FF(my_status)
582	–	if (my_status .ne. 0) then
583	–	thisStat = -1
584	–	haveSaneForceField = .false.
585	–	return
586	–	endif
587	–	endif
588	–
764		if (.not. haveNeighborList) then
765		!! Create neighbor lists
766		call expandNeighborList(nLocal, my_status)
#	Line 619 | Line 794 | contains
794
795		!! Stress Tensor
796		real( kind = dp), dimension(9) :: tau
797	<	real ( kind = dp ) :: pot
797	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
798		logical ( kind = 2) :: do_pot_c, do_stress_c
799		logical :: do_pot
800		logical :: do_stress
801		logical :: in_switching_region
802		#ifdef IS_MPI
803	<	real( kind = DP ) :: pot_local
803	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
804		integer :: nAtomsInRow
805		integer :: nAtomsInCol
806		integer :: nprocs
#	Line 640 | Line 815 | contains
815		integer :: nlist
816		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
817		real( kind = DP ) :: sw, dswdr, swderiv, mf
818	+	real( kind = DP ) :: rVal
819		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
820		real(kind=dp) :: rfpot, mu_i, virial
821	+	real(kind=dp):: rCut
822		integer :: me_i, me_j, n_in_i, n_in_j
823		logical :: is_dp_i
824		integer :: neighborListSize
#	Line 650 | Line 827 | contains
827		integer :: propPack_i, propPack_j
828		integer :: loopStart, loopEnd, loop
829		integer :: iHash
830	<	real(kind=dp) :: listSkin = 1.0
830	>	integer :: i1
831
832	+	!! the variables for the box dipole moment
833	+	#ifdef IS_MPI
834	+	integer :: pChgCount_local
835	+	integer :: nChgCount_local
836	+	real(kind=dp) :: pChg_local
837	+	real(kind=dp) :: nChg_local
838	+	real(kind=dp), dimension(3) :: pChgPos_local
839	+	real(kind=dp), dimension(3) :: nChgPos_local
840	+	real(kind=dp), dimension(3) :: dipVec_local
841	+	#endif
842	+	integer :: pChgCount
843	+	integer :: nChgCount
844	+	real(kind=dp) :: pChg
845	+	real(kind=dp) :: nChg
846	+	real(kind=dp) :: chg_value
847	+	real(kind=dp), dimension(3) :: pChgPos
848	+	real(kind=dp), dimension(3) :: nChgPos
849	+	real(kind=dp), dimension(3) :: dipVec
850	+	real(kind=dp), dimension(3) :: chgVec
851	+
852	+	!! initialize box dipole variables
853	+	if (do_box_dipole) then
854	+	#ifdef IS_MPI
855	+	pChg_local = 0.0_dp
856	+	nChg_local = 0.0_dp
857	+	pChgCount_local = 0
858	+	nChgCount_local = 0
859	+	do i=1, 3
860	+	pChgPos_local = 0.0_dp
861	+	nChgPos_local = 0.0_dp
862	+	dipVec_local = 0.0_dp
863	+	enddo
864	+	#endif
865	+	pChg = 0.0_dp
866	+	nChg = 0.0_dp
867	+	pChgCount = 0
868	+	nChgCount = 0
869	+	chg_value = 0.0_dp
870	+
871	+	do i=1, 3
872	+	pChgPos(i) = 0.0_dp
873	+	nChgPos(i) = 0.0_dp
874	+	dipVec(i) = 0.0_dp
875	+	chgVec(i) = 0.0_dp
876	+	boxDipole(i) = 0.0_dp
877	+	enddo
878	+	endif
879	+
880		!! initialize local variables
881
882		#ifdef IS_MPI
#	Line 714 | Line 939 | contains
939		! (but only on the first time through):
940		if (loop .eq. loopStart) then
941		#ifdef IS_MPI
942	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
942	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
943		update_nlist)
944		#else
945	<	call checkNeighborList(nGroups, q_group, listSkin, &
945	>	call checkNeighborList(nGroups, q_group, skinThickness, &
946		update_nlist)
947		#endif
948		endif
#	Line 775 | Line 1000 | contains
1000		me_j = atid(j)
1001		call get_interatomic_vector(q_group(:,i), &
1002		q_group(:,j), d_grp, rgrpsq)
1003	<	#endif
1003	>	#endif
1004
1005	<	if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
1005	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
1006		if (update_nlist) then
1007		nlist = nlist + 1
1008
#	Line 797 | Line 1022 | contains
1022
1023		list(nlist) = j
1024		endif
1025	+
1026	+	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
1027
1028	<	if (loop .eq. PAIR_LOOP) then
1029	<	vij = 0.0d0
1030	<	fij(1:3) = 0.0d0
1031	<	endif
1032	<
1033	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
1034	<	in_switching_region)
1035	<
1036	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
1037	<
1038	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
1039	<
1040	<	atom1 = groupListRow(ia)
1041	<
1042	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1043	<
1044	<	atom2 = groupListCol(jb)
1045	<
1046	<	if (skipThisPair(atom1, atom2)) cycle inner
1047	<
1048	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1049	<	d_atm(1:3) = d_grp(1:3)
1050	<	ratmsq = rgrpsq
1051	<	else
1028	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
1029	>	if (loop .eq. PAIR_LOOP) then
1030	>	vij = 0.0_dp
1031	>	fij(1) = 0.0_dp
1032	>	fij(2) = 0.0_dp
1033	>	fij(3) = 0.0_dp
1034	>	endif
1035	>
1036	>	call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
1037	>
1038	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
1039	>
1040	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
1041	>
1042	>	atom1 = groupListRow(ia)
1043	>
1044	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1045	>
1046	>	atom2 = groupListCol(jb)
1047	>
1048	>	if (skipThisPair(atom1, atom2))  cycle inner
1049	>
1050	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1051	>	d_atm(1) = d_grp(1)
1052	>	d_atm(2) = d_grp(2)
1053	>	d_atm(3) = d_grp(3)
1054	>	ratmsq = rgrpsq
1055	>	else
1056		#ifdef IS_MPI
1057	<	call get_interatomic_vector(q_Row(:,atom1), &
1058	<	q_Col(:,atom2), d_atm, ratmsq)
1057	>	call get_interatomic_vector(q_Row(:,atom1), &
1058	>	q_Col(:,atom2), d_atm, ratmsq)
1059		#else
1060	<	call get_interatomic_vector(q(:,atom1), &
1061	<	q(:,atom2), d_atm, ratmsq)
1060	>	call get_interatomic_vector(q(:,atom1), &
1061	>	q(:,atom2), d_atm, ratmsq)
1062		#endif
1063	<	endif
1064	<
1065	<	if (loop .eq. PREPAIR_LOOP) then
1063	>	endif
1064	>
1065	>	if (loop .eq. PREPAIR_LOOP) then
1066		#ifdef IS_MPI
1067	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1068	<	rgrpsq, d_grp, do_pot, do_stress, &
1069	<	eFrame, A, f, t, pot_local)
1067	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1068	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1069	>	eFrame, A, f, t, pot_local)
1070		#else
1071	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1072	<	rgrpsq, d_grp, do_pot, do_stress, &
1073	<	eFrame, A, f, t, pot)
1071	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1072	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1073	>	eFrame, A, f, t, pot)
1074		#endif
1075	<	else
1075	>	else
1076		#ifdef IS_MPI
1077	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1078	<	do_pot, &
1079	<	eFrame, A, f, t, pot_local, vpair, fpair)
1077	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1078	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1079	>	fpair, d_grp, rgrp, rCut)
1080		#else
1081	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1082	<	do_pot,  &
1083	<	eFrame, A, f, t, pot, vpair, fpair)
1081	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1082	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1083	>	d_grp, rgrp, rCut)
1084		#endif
1085	+	vij = vij + vpair
1086	+	fij(1) = fij(1) + fpair(1)
1087	+	fij(2) = fij(2) + fpair(2)
1088	+	fij(3) = fij(3) + fpair(3)
1089	+	endif
1090	+	enddo inner
1091	+	enddo
1092
1093	<	vij = vij + vpair
1094	<	fij(1:3) = fij(1:3) + fpair(1:3)
1095	<	endif
1096	<	enddo inner
1097	<	enddo
1098	<
1099	<	if (loop .eq. PAIR_LOOP) then
1100	<	if (in_switching_region) then
1101	<	swderiv = vij*dswdr/rgrp
1102	<	fij(1) = fij(1) + swderiv*d_grp(1)
865	<	fij(2) = fij(2) + swderiv*d_grp(2)
866	<	fij(3) = fij(3) + swderiv*d_grp(3)
867	<
868	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
869	<	atom1=groupListRow(ia)
870	<	mf = mfactRow(atom1)
1093	>	if (loop .eq. PAIR_LOOP) then
1094	>	if (in_switching_region) then
1095	>	swderiv = vij*dswdr/rgrp
1096	>	fij(1) = fij(1) + swderiv*d_grp(1)
1097	>	fij(2) = fij(2) + swderiv*d_grp(2)
1098	>	fij(3) = fij(3) + swderiv*d_grp(3)
1099	>
1100	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1101	>	atom1=groupListRow(ia)
1102	>	mf = mfactRow(atom1)
1103		#ifdef IS_MPI
1104	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1105	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1106	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1104	>	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1105	>	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1106	>	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1107		#else
1108	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1109	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1110	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1108	>	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1109	>	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1110	>	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1111		#endif
1112	<	enddo
1113	<
1114	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1115	<	atom2=groupListCol(jb)
1116	<	mf = mfactCol(atom2)
1112	>	enddo
1113	>
1114	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1115	>	atom2=groupListCol(jb)
1116	>	mf = mfactCol(atom2)
1117		#ifdef IS_MPI
1118	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1119	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1120	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1118	>	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1119	>	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1120	>	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1121		#else
1122	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1123	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1124	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1122	>	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1123	>	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1124	>	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1125		#endif
1126	<	enddo
1127	<	endif
1126	>	enddo
1127	>	endif
1128
1129	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1129	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1130	>	endif
1131		endif
1132	<	end if
1132	>	endif
1133		enddo
1134	+
1135		enddo outer
1136
1137		if (update_nlist) then
#	Line 957 | Line 1191 | contains
1191
1192		if (do_pot) then
1193		! scatter/gather pot_row into the members of my column
1194	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1195	<
1194	>	do i = 1,LR_POT_TYPES
1195	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1196	>	end do
1197		! scatter/gather pot_local into all other procs
1198		! add resultant to get total pot
1199		do i = 1, nlocal
1200	<	pot_local = pot_local + pot_Temp(i)
1200	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1201	>	+ pot_Temp(1:LR_POT_TYPES,i)
1202		enddo
1203
1204		pot_Temp = 0.0_DP
1205	<
1206	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1205	>	do i = 1,LR_POT_TYPES
1206	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1207	>	end do
1208		do i = 1, nlocal
1209	<	pot_local = pot_local + pot_Temp(i)
1209	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1210	>	+ pot_Temp(1:LR_POT_TYPES,i)
1211		enddo
1212
1213		endif
1214		#endif
1215
1216	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1216	>	if (SIM_requires_postpair_calc) then
1217	>	do i = 1, nlocal
1218	>
1219	>	! we loop only over the local atoms, so we don't need row and column
1220	>	! lookups for the types
1221	>
1222	>	me_i = atid(i)
1223	>
1224	>	! is the atom electrostatic?  See if it would have an
1225	>	! electrostatic interaction with itself
1226	>	iHash = InteractionHash(me_i,me_i)
1227
1228	<	if (FF_uses_RF .and. SIM_uses_RF) then
981	<
1228	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1229		#ifdef IS_MPI
1230	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1231	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1232	<	do i = 1,nlocal
1233	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1234	<	end do
1230	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1231	>	t, do_pot)
1232	>	#else
1233	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1234	>	t, do_pot)
1235		#endif
1236	<
1237	<	do i = 1, nLocal
1238	<
1239	<	rfpot = 0.0_DP
1236	>	endif
1237	>
1238	>
1239	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1240	>
1241	>	! loop over the excludes to accumulate RF stuff we've
1242	>	! left out of the normal pair loop
1243	>
1244	>	do i1 = 1, nSkipsForAtom(i)
1245	>	j = skipsForAtom(i, i1)
1246	>
1247	>	! prevent overcounting of the skips
1248	>	if (i.lt.j) then
1249	>	call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1250	>	rVal = sqrt(ratmsq)
1251	>	call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1252		#ifdef IS_MPI
1253	<	me_i = atid_row(i)
1253	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1254	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1255		#else
1256	<	me_i = atid(i)
1256	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1257	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1258		#endif
1259	<	iHash = InteractionHash(me_i,me_j)
1260	<
1261	<	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1259	>	endif
1260	>	enddo
1261	>	endif
1262
1263	<	mu_i = getDipoleMoment(me_i)
1003	<
1004	<	!! The reaction field needs to include a self contribution
1005	<	!! to the field:
1006	<	call accumulate_self_rf(i, mu_i, eFrame)
1007	<	!! Get the reaction field contribution to the
1008	<	!! potential and torques:
1009	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1263	>	if (do_box_dipole) then
1264		#ifdef IS_MPI
1265	<	pot_local = pot_local + rfpot
1265	>	call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1266	>	nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1267	>	pChgCount_local, nChgCount_local)
1268		#else
1269	<	pot = pot + rfpot
1270	<
1269	>	call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1270	>	pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1271		#endif
1272	<	endif
1273	<	enddo
1018	<	endif
1272	>	endif
1273	>	enddo
1274		endif
1275
1021	–
1276		#ifdef IS_MPI
1023	–
1277		if (do_pot) then
1278	<	pot = pot + pot_local
1279	<	!! we assume the c code will do the allreduce to get the total potential
1280	<	!! we could do it right here if we needed to...
1278	>	#ifdef SINGLE_PRECISION
1279	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1280	>	mpi_comm_world,mpi_err)
1281	>	#else
1282	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1283	>	mpi_sum, mpi_comm_world,mpi_err)
1284	>	#endif
1285		endif
1286	<
1286	>
1287		if (do_stress) then
1288	+	#ifdef SINGLE_PRECISION
1289	+	call mpi_allreduce(tau_Temp, tau, 9,mpi_real,mpi_sum, &
1290	+	mpi_comm_world,mpi_err)
1291	+	call mpi_allreduce(virial_Temp, virial,1,mpi_real,mpi_sum, &
1292	+	mpi_comm_world,mpi_err)
1293	+	#else
1294		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1295		mpi_comm_world,mpi_err)
1296		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1297		mpi_comm_world,mpi_err)
1298	+	#endif
1299		endif
1300	+
1301	+	if (do_box_dipole) then
1302
1303	+	#ifdef SINGLE_PRECISION
1304	+	call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1305	+	mpi_comm_world, mpi_err)
1306	+	call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1307	+	mpi_comm_world, mpi_err)
1308	+	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1309	+	mpi_comm_world, mpi_err)
1310	+	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1311	+	mpi_comm_world, mpi_err)
1312	+	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1313	+	mpi_comm_world, mpi_err)
1314	+	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1315	+	mpi_comm_world, mpi_err)
1316	+	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1317	+	mpi_comm_world, mpi_err)
1318		#else
1319	+	call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1320	+	mpi_comm_world, mpi_err)
1321	+	call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1322	+	mpi_comm_world, mpi_err)
1323	+	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1324	+	mpi_sum, mpi_comm_world, mpi_err)
1325	+	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1326	+	mpi_sum, mpi_comm_world, mpi_err)
1327	+	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1328	+	mpi_sum, mpi_comm_world, mpi_err)
1329	+	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1330	+	mpi_sum, mpi_comm_world, mpi_err)
1331	+	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1332	+	mpi_sum, mpi_comm_world, mpi_err)
1333	+	#endif
1334
1335	+	endif
1336	+
1337	+	#else
1338	+
1339		if (do_stress) then
1340		tau = tau_Temp
1341		virial = virial_Temp
1342		endif
1343	<
1343	>
1344		#endif
1345
1346	+	if (do_box_dipole) then
1347	+	! first load the accumulated dipole moment (if dipoles were present)
1348	+	boxDipole(1) = dipVec(1)
1349	+	boxDipole(2) = dipVec(2)
1350	+	boxDipole(3) = dipVec(3)
1351	+
1352	+	! now include the dipole moment due to charges
1353	+	! use the lesser of the positive and negative charge totals
1354	+	if (nChg .le. pChg) then
1355	+	chg_value = nChg
1356	+	else
1357	+	chg_value = pChg
1358	+	endif
1359	+
1360	+	! find the average positions
1361	+	if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1362	+	pChgPos = pChgPos / pChgCount
1363	+	nChgPos = nChgPos / nChgCount
1364	+	endif
1365	+
1366	+	! dipole is from the negative to the positive (physics notation)
1367	+	chgVec(1) = pChgPos(1) - nChgPos(1)
1368	+	chgVec(2) = pChgPos(2) - nChgPos(2)
1369	+	chgVec(3) = pChgPos(3) - nChgPos(3)
1370	+
1371	+	boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1372	+	boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1373	+	boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1374	+
1375	+	endif
1376	+
1377		end subroutine do_force_loop
1378
1379		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1380	<	eFrame, A, f, t, pot, vpair, fpair)
1380	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1381
1382	<	real( kind = dp ) :: pot, vpair, sw
1382	>	real( kind = dp ) :: vpair, sw
1383	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1384		real( kind = dp ), dimension(3) :: fpair
1385		real( kind = dp ), dimension(nLocal)   :: mfact
1386		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1059 | Line 1391 | contains
1391		logical, intent(inout) :: do_pot
1392		integer, intent(in) :: i, j
1393		real ( kind = dp ), intent(inout) :: rijsq
1394	<	real ( kind = dp )                :: r
1394	>	real ( kind = dp ), intent(inout) :: r_grp
1395		real ( kind = dp ), intent(inout) :: d(3)
1396	<	real ( kind = dp ) :: ebalance
1396	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1397	>	real ( kind = dp ), intent(inout) :: rCut
1398	>	real ( kind = dp ) :: r
1399	>	real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1400		integer :: me_i, me_j
1401	+	integer :: k
1402
1403		integer :: iHash
1404
1405		r = sqrt(rijsq)
1406	<	vpair = 0.0d0
1407	<	fpair(1:3) = 0.0d0
1406	>
1407	>	vpair = 0.0_dp
1408	>	fpair(1:3) = 0.0_dp
1409
1410		#ifdef IS_MPI
1411		me_i = atid_row(i)
#	Line 1079 | Line 1416 | contains
1416		#endif
1417
1418		iHash = InteractionHash(me_i, me_j)
1419	<
1419	>
1420		if ( iand(iHash, LJ_PAIR).ne.0 ) then
1421	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1421	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1422	>	pot(VDW_POT), f, do_pot)
1423		endif
1424	<
1424	>
1425		if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1426	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1427	<	pot, eFrame, f, t, do_pot, corrMethod)
1090	<
1091	<	if (FF_uses_RF .and. SIM_uses_RF) then
1092	<
1093	<	! CHECK ME (RF needs to know about all electrostatic types)
1094	<	call accumulate_rf(i, j, r, eFrame, sw)
1095	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1096	<	endif
1097	<
1426	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1427	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1428		endif
1429	<
1429	>
1430		if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1431		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1432	<	pot, A, f, t, do_pot)
1432	>	pot(HB_POT), A, f, t, do_pot)
1433		endif
1434	<
1434	>
1435		if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1436		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1437	<	pot, A, f, t, do_pot)
1437	>	pot(HB_POT), A, f, t, do_pot)
1438		endif
1439	<
1439	>
1440		if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1441		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1442	<	pot, A, f, t, do_pot)
1442	>	pot(VDW_POT), A, f, t, do_pot)
1443		endif
1444
1445		if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1446	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1447	<	!           pot, A, f, t, do_pot)
1446	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1447	>	pot(VDW_POT), A, f, t, do_pot)
1448		endif
1449	<
1449	>
1450		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1451	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1452	<	do_pot)
1451	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1452	>	pot(METALLIC_POT), f, do_pot)
1453		endif
1454	<
1454	>
1455		if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1456		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1457	<	pot, A, f, t, do_pot)
1457	>	pot(VDW_POT), A, f, t, do_pot)
1458		endif
1459	<
1459	>
1460		if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1461		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1462	<	pot, A, f, t, do_pot)
1462	>	pot(VDW_POT), A, f, t, do_pot)
1463		endif
1464	<
1464	>
1465	>	if ( iand(iHash, SC_PAIR).ne.0 ) then
1466	>	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1467	>	pot(METALLIC_POT), f, do_pot)
1468	>	endif
1469	>
1470		end subroutine do_pair
1471
1472	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1472	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1473		do_pot, do_stress, eFrame, A, f, t, pot)
1474
1475	<	real( kind = dp ) :: pot, sw
1475	>	real( kind = dp ) :: sw
1476	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1477		real( kind = dp ), dimension(9,nLocal) :: eFrame
1478		real (kind=dp), dimension(9,nLocal) :: A
1479		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1145 | Line 1481 | contains
1481
1482		logical, intent(inout) :: do_pot, do_stress
1483		integer, intent(in) :: i, j
1484	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1484	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1485		real ( kind = dp )                :: r, rc
1486		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1487
1488		integer :: me_i, me_j, iHash
1489
1490	+	r = sqrt(rijsq)
1491	+
1492		#ifdef IS_MPI
1493		me_i = atid_row(i)
1494		me_j = atid_col(j)
#	Line 1162 | Line 1500 | contains
1500		iHash = InteractionHash(me_i, me_j)
1501
1502		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1503	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1503	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1504		endif
1505	+
1506	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1507	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1508	+	endif
1509
1510		end subroutine do_prepair
1511
1512
1513		subroutine do_preforce(nlocal,pot)
1514		integer :: nlocal
1515	<	real( kind = dp ) :: pot
1515	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1516
1517		if (FF_uses_EAM .and. SIM_uses_EAM) then
1518	<	call calc_EAM_preforce_Frho(nlocal,pot)
1518	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1519		endif
1520	<
1521	<
1520	>	if (FF_uses_SC .and. SIM_uses_SC) then
1521	>	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1522	>	endif
1523		end subroutine do_preforce
1524
1525
#	Line 1188 | Line 1531 | contains
1531		real( kind = dp ) :: d(3), scaled(3)
1532		integer i
1533
1534	<	d(1:3) = q_j(1:3) - q_i(1:3)
1534	>	d(1) = q_j(1) - q_i(1)
1535	>	d(2) = q_j(2) - q_i(2)
1536	>	d(3) = q_j(3) - q_i(3)
1537
1538		! Wrap back into periodic box if necessary
1539		if ( SIM_uses_PBC ) then
1540
1541		if( .not.boxIsOrthorhombic ) then
1542		! calc the scaled coordinates.
1543	+	! scaled = matmul(HmatInv, d)
1544
1545	<	scaled = matmul(HmatInv, d)
1546	<
1545	>	scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1546	>	scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1547	>	scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1548	>
1549		! wrap the scaled coordinates
1550
1551	<	scaled = scaled  - anint(scaled)
1551	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1552	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1553	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1554
1205	–
1555		! calc the wrapped real coordinates from the wrapped scaled
1556		! coordinates
1557	+	! d = matmul(Hmat,scaled)
1558	+	d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1559	+	d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1560	+	d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1561
1209	–	d = matmul(Hmat,scaled)
1210	–
1562		else
1563		! calc the scaled coordinates.
1564
1565	<	do i = 1, 3
1566	<	scaled(i) = d(i) * HmatInv(i,i)
1565	>	scaled(1) = d(1) * HmatInv(1,1)
1566	>	scaled(2) = d(2) * HmatInv(2,2)
1567	>	scaled(3) = d(3) * HmatInv(3,3)
1568	>
1569	>	! wrap the scaled coordinates
1570	>
1571	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1572	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1573	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1574
1575	<	! wrap the scaled coordinates
1575	>	! calc the wrapped real coordinates from the wrapped scaled
1576	>	! coordinates
1577
1578	<	scaled(i) = scaled(i) - anint(scaled(i))
1578	>	d(1) = scaled(1)*Hmat(1,1)
1579	>	d(2) = scaled(2)*Hmat(2,2)
1580	>	d(3) = scaled(3)*Hmat(3,3)
1581
1221	–	! calc the wrapped real coordinates from the wrapped scaled
1222	–	! coordinates
1223	–
1224	–	d(i) = scaled(i)*Hmat(i,i)
1225	–	enddo
1582		endif
1583
1584		endif
1585
1586	<	r_sq = dot_product(d,d)
1586	>	r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1587
1588		end subroutine get_interatomic_vector
1589
#	Line 1259 | Line 1615 | contains
1615		pot_Col = 0.0_dp
1616		pot_Temp = 0.0_dp
1617
1262	–	rf_Row = 0.0_dp
1263	–	rf_Col = 0.0_dp
1264	–	rf_Temp = 0.0_dp
1265	–
1618		#endif
1619
1620		if (FF_uses_EAM .and. SIM_uses_EAM) then
1621		call clean_EAM()
1622		endif
1623
1272	–	rf = 0.0_dp
1624		tau_Temp = 0.0_dp
1625		virial_Temp = 0.0_dp
1626		end subroutine zero_work_arrays
#	Line 1363 | Line 1714 | contains
1714
1715		function FF_RequiresPrepairCalc() result(doesit)
1716		logical :: doesit
1717	<	doesit = FF_uses_EAM
1717	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1718	>	.or. FF_uses_MEAM
1719		end function FF_RequiresPrepairCalc
1720
1369	–	function FF_RequiresPostpairCalc() result(doesit)
1370	–	logical :: doesit
1371	–	doesit = FF_uses_RF
1372	–	end function FF_RequiresPostpairCalc
1373	–
1721		#ifdef PROFILE
1722		function getforcetime() result(totalforcetime)
1723		real(kind=dp) :: totalforcetime

Diff Legend

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