ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE-4/src/UseTheForce/doForces.F90

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2309 by chrisfen, Sun Sep 18 20:45:38 2005 UTC vs.
Revision 3133 by chuckv, Tue May 22 19:30:27 2007 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.46 2005-09-18 20:45:38 chrisfen Exp $, $Date: 2005-09-18 20:45:38 $, $Name: not supported by cvs2svn $, $Revision: 1.46 $
48	>	!! @version $Id: doForces.F90,v 1.90 2007-05-22 19:30:27 chuckv Exp $, $Date: 2007-05-22 19:30:27 $, $Name: not supported by cvs2svn $, $Revision: 1.90 $
49
50
51		module doForces
#	Line 58 | Line 58 | module doForces
58		use lj
59		use sticky
60		use electrostatic_module
61	<	use reaction_field_module
62	<	use gb_pair
61	>	use gayberne
62		use shapes
63		use vector_class
64		use eam
65	+	use suttonchen
66		use status
67		#ifdef IS_MPI
68		use mpiSimulation
#	Line 72 | Line 72 | module doForces
72		PRIVATE
73
74		#define __FORTRAN90
75	–	#include "UseTheForce/fSwitchingFunction.h"
75		#include "UseTheForce/fCutoffPolicy.h"
76		#include "UseTheForce/DarkSide/fInteractionMap.h"
77	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78
79	–
79		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
80		INTEGER, PARAMETER:: PAIR_LOOP    = 2
81
#	Line 85 | Line 84 | module doForces
84		logical, save :: haveSaneForceField = .false.
85		logical, save :: haveInteractionHash = .false.
86		logical, save :: haveGtypeCutoffMap = .false.
87	<	logical, save :: haveRlist = .false.
87	>	logical, save :: haveDefaultCutoffs = .false.
88	>	logical, save :: haveSkinThickness = .false.
89	>	logical, save :: haveElectrostaticSummationMethod = .false.
90	>	logical, save :: haveCutoffPolicy = .false.
91	>	logical, save :: VisitCutoffsAfterComputing = .false.
92	>	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_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_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	+	logical, save :: SIM_uses_AtomicVirial
110
111		integer, save :: electrostaticSummationMethod
112	+	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
113
114	+	real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
115	+	real(kind=dp), save :: skinThickness
116	+	logical, save :: defaultDoShiftPot
117	+	logical, save :: defaultDoShiftFrc
118	+
119		public :: init_FF
120	<	public :: setDefaultCutoffs
120	>	public :: setCutoffs
121	>	public :: cWasLame
122	>	public :: setElectrostaticMethod
123	>	public :: setBoxDipole
124	>	public :: getBoxDipole
125	>	public :: setCutoffPolicy
126	>	public :: setSkinThickness
127		public :: do_force_loop
106	–	public :: createInteractionHash
107	–	public :: createGtypeCutoffMap
108	–	public :: getStickyCut
109	–	public :: getStickyPowerCut
110	–	public :: getGayBerneCut
111	–	public :: getEAMCut
112	–	public :: getShapeCut
128
129		#ifdef PROFILE
130		public :: getforcetime
#	Line 122 | Line 137 | module doForces
137		! Bit hash to determine pair-pair interactions.
138		integer, dimension(:,:), allocatable :: InteractionHash
139		real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
140	<	real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
141	<	integer, dimension(:), allocatable :: groupToGtype
142	<	real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
140	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
141	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
142	>
143	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
144	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
145	>
146	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
147	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
148		type ::gtypeCutoffs
149		real(kind=dp) :: rcut
150		real(kind=dp) :: rcutsq
#	Line 132 | Line 152 | module doForces
152		end type gtypeCutoffs
153		type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
154
155	<	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
156	<	real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
137	<
155	>	real(kind=dp), dimension(3) :: boxDipole
156	>
157		contains
158
159	<	subroutine createInteractionHash(status)
159	>	subroutine createInteractionHash()
160		integer :: nAtypes
142	–	integer, intent(out) :: status
161		integer :: i
162		integer :: j
163		integer :: iHash
#	Line 151 | Line 169 | contains
169		logical :: i_is_GB
170		logical :: i_is_EAM
171		logical :: i_is_Shape
172	+	logical :: i_is_SC
173	+	logical :: i_is_MEAM
174		logical :: j_is_LJ
175		logical :: j_is_Elect
176		logical :: j_is_Sticky
#	Line 158 | Line 178 | contains
178		logical :: j_is_GB
179		logical :: j_is_EAM
180		logical :: j_is_Shape
181	+	logical :: j_is_SC
182	+	logical :: j_is_MEAM
183		real(kind=dp) :: myRcut
184
163	–	status = 0
164	–
185		if (.not. associated(atypes)) then
186	<	call handleError("atype", "atypes was not present before call of createInteractionHash!")
167	<	status = -1
186	>	call handleError("doForces", "atypes was not present before call of createInteractionHash!")
187		return
188		endif
189
190		nAtypes = getSize(atypes)
191
192		if (nAtypes == 0) then
193	<	status = -1
193	>	call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
194		return
195		end if
196
197		if (.not. allocated(InteractionHash)) then
198		allocate(InteractionHash(nAtypes,nAtypes))
199	+	else
200	+	deallocate(InteractionHash)
201	+	allocate(InteractionHash(nAtypes,nAtypes))
202		endif
203
204		if (.not. allocated(atypeMaxCutoff)) then
205		allocate(atypeMaxCutoff(nAtypes))
206	+	else
207	+	deallocate(atypeMaxCutoff)
208	+	allocate(atypeMaxCutoff(nAtypes))
209		endif
210
211		do i = 1, nAtypes
#	Line 191 | Line 216 | contains
216		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
217		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
218		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
219	+	call getElementProperty(atypes, i, "is_SC", i_is_SC)
220	+	call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
221
222		do j = i, nAtypes
223
#	Line 204 | Line 231 | contains
231		call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
232		call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
233		call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
234	+	call getElementProperty(atypes, j, "is_SC", j_is_SC)
235	+	call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
236
237		if (i_is_LJ .and. j_is_LJ) then
238		iHash = ior(iHash, LJ_PAIR)
#	Line 223 | Line 252 | contains
252
253		if (i_is_EAM .and. j_is_EAM) then
254		iHash = ior(iHash, EAM_PAIR)
255	+	endif
256	+
257	+	if (i_is_SC .and. j_is_SC) then
258	+	iHash = ior(iHash, SC_PAIR)
259		endif
260
261		if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
#	Line 244 | Line 277 | contains
277		haveInteractionHash = .true.
278		end subroutine createInteractionHash
279
280	<	subroutine createGtypeCutoffMap(stat)
280	>	subroutine createGtypeCutoffMap()
281
249	–	integer, intent(out), optional :: stat
282		logical :: i_is_LJ
283		logical :: i_is_Elect
284		logical :: i_is_Sticky
#	Line 254 | Line 286 | contains
286		logical :: i_is_GB
287		logical :: i_is_EAM
288		logical :: i_is_Shape
289	+	logical :: i_is_SC
290		logical :: GtypeFound
291
292		integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
293	<	integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
293	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
294		integer :: nGroupsInRow
295	<	real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
295	>	integer :: nGroupsInCol
296	>	integer :: nGroupTypesRow,nGroupTypesCol
297	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
298		real(kind=dp) :: biggestAtypeCutoff
299
265	–	stat = 0
300		if (.not. haveInteractionHash) then
301	<	call createInteractionHash(myStatus)
268	<	if (myStatus .ne. 0) then
269	<	write(default_error, *) 'createInteractionHash failed in doForces!'
270	<	stat = -1
271	<	return
272	<	endif
301	>	call createInteractionHash()
302		endif
303		#ifdef IS_MPI
304		nGroupsInRow = getNgroupsInRow(plan_group_row)
305	+	nGroupsInCol = getNgroupsInCol(plan_group_col)
306		#endif
307		nAtypes = getSize(atypes)
308		! Set all of the initial cutoffs to zero.
#	Line 286 | Line 316 | contains
316		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
317		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
318		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
319	<
319	>	call getElementProperty(atypes, i, "is_SC", i_is_SC)
320
321	<	if (i_is_LJ) then
322	<	thisRcut = getSigma(i) * 2.5_dp
323	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
324	<	endif
325	<	if (i_is_Elect) then
326	<	thisRcut = defaultRcut
327	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
328	<	endif
329	<	if (i_is_Sticky) then
330	<	thisRcut = getStickyCut(i)
331	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
332	<	endif
333	<	if (i_is_StickyP) then
334	<	thisRcut = getStickyPowerCut(i)
335	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
336	<	endif
337	<	if (i_is_GB) then
338	<	thisRcut = getGayBerneCut(i)
339	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
340	<	endif
341	<	if (i_is_EAM) then
342	<	thisRcut = getEAMCut(i)
343	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
344	<	endif
345	<	if (i_is_Shape) then
346	<	thisRcut = getShapeCut(i)
347	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321	>	if (haveDefaultCutoffs) then
322	>	atypeMaxCutoff(i) = defaultRcut
323	>	else
324	>	if (i_is_LJ) then
325	>	thisRcut = getSigma(i) * 2.5_dp
326	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
327	>	endif
328	>	if (i_is_Elect) then
329	>	thisRcut = defaultRcut
330	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
331	>	endif
332	>	if (i_is_Sticky) then
333	>	thisRcut = getStickyCut(i)
334	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
335	>	endif
336	>	if (i_is_StickyP) then
337	>	thisRcut = getStickyPowerCut(i)
338	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
339	>	endif
340	>	if (i_is_GB) then
341	>	thisRcut = getGayBerneCut(i)
342	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
343	>	endif
344	>	if (i_is_EAM) then
345	>	thisRcut = getEAMCut(i)
346	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
347	>	endif
348	>	if (i_is_Shape) then
349	>	thisRcut = getShapeCut(i)
350	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
351	>	endif
352	>	if (i_is_SC) then
353	>	thisRcut = getSCCut(i)
354	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
355	>	endif
356		endif
357	<
357	>
358		if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
359		biggestAtypeCutoff = atypeMaxCutoff(i)
360		endif
361	+
362		endif
363		enddo
325	–
326	–	nGroupTypes = 0
364
365		istart = 1
366	+	jstart = 1
367		#ifdef IS_MPI
368		iend = nGroupsInRow
369	+	jend = nGroupsInCol
370		#else
371		iend = nGroups
372	+	jend = nGroups
373		#endif
374
375		!! allocate the groupToGtype and gtypeMaxCutoff here.
376	<	if(.not.allocated(groupToGtype)) then
377	<	allocate(groupToGtype(iend))
378	<	allocate(groupMaxCutoff(iend))
379	<	allocate(gtypeMaxCutoff(iend))
380	<	groupMaxCutoff = 0.0_dp
381	<	gtypeMaxCutoff = 0.0_dp
376	>	if(.not.allocated(groupToGtypeRow)) then
377	>	!  allocate(groupToGtype(iend))
378	>	allocate(groupToGtypeRow(iend))
379	>	else
380	>	deallocate(groupToGtypeRow)
381	>	allocate(groupToGtypeRow(iend))
382		endif
383	+	if(.not.allocated(groupMaxCutoffRow)) then
384	+	allocate(groupMaxCutoffRow(iend))
385	+	else
386	+	deallocate(groupMaxCutoffRow)
387	+	allocate(groupMaxCutoffRow(iend))
388	+	end if
389	+
390	+	if(.not.allocated(gtypeMaxCutoffRow)) then
391	+	allocate(gtypeMaxCutoffRow(iend))
392	+	else
393	+	deallocate(gtypeMaxCutoffRow)
394	+	allocate(gtypeMaxCutoffRow(iend))
395	+	endif
396	+
397	+
398	+	#ifdef IS_MPI
399	+	! We only allocate new storage if we are in MPI because Ncol /= Nrow
400	+	if(.not.associated(groupToGtypeCol)) then
401	+	allocate(groupToGtypeCol(jend))
402	+	else
403	+	deallocate(groupToGtypeCol)
404	+	allocate(groupToGtypeCol(jend))
405	+	end if
406	+
407	+	if(.not.associated(groupMaxCutoffCol)) then
408	+	allocate(groupMaxCutoffCol(jend))
409	+	else
410	+	deallocate(groupMaxCutoffCol)
411	+	allocate(groupMaxCutoffCol(jend))
412	+	end if
413	+	if(.not.associated(gtypeMaxCutoffCol)) then
414	+	allocate(gtypeMaxCutoffCol(jend))
415	+	else
416	+	deallocate(gtypeMaxCutoffCol)
417	+	allocate(gtypeMaxCutoffCol(jend))
418	+	end if
419	+
420	+	groupMaxCutoffCol = 0.0_dp
421	+	gtypeMaxCutoffCol = 0.0_dp
422	+
423	+	#endif
424	+	groupMaxCutoffRow = 0.0_dp
425	+	gtypeMaxCutoffRow = 0.0_dp
426	+
427	+
428		!! first we do a single loop over the cutoff groups to find the
429		!! largest cutoff for any atypes present in this group.  We also
430		!! create gtypes at this point.
431
432	<	tol = 1.0d-6
433	<
432	>	tol = 1.0e-6_dp
433	>	nGroupTypesRow = 0
434	>	nGroupTypesCol = 0
435		do i = istart, iend
436		n_in_i = groupStartRow(i+1) - groupStartRow(i)
437	<	groupMaxCutoff(i) = 0.0_dp
437	>	groupMaxCutoffRow(i) = 0.0_dp
438		do ia = groupStartRow(i), groupStartRow(i+1)-1
439		atom1 = groupListRow(ia)
440		#ifdef IS_MPI
#	Line 356 | Line 442 | contains
442		#else
443		me_i = atid(atom1)
444		#endif
445	<	if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
446	<	groupMaxCutoff(i)=atypeMaxCutoff(me_i)
445	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
446	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
447		endif
448		enddo
449	+	if (nGroupTypesRow.eq.0) then
450	+	nGroupTypesRow = nGroupTypesRow + 1
451	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
452	+	groupToGtypeRow(i) = nGroupTypesRow
453	+	else
454	+	GtypeFound = .false.
455	+	do g = 1, nGroupTypesRow
456	+	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
457	+	groupToGtypeRow(i) = g
458	+	GtypeFound = .true.
459	+	endif
460	+	enddo
461	+	if (.not.GtypeFound) then
462	+	nGroupTypesRow = nGroupTypesRow + 1
463	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
464	+	groupToGtypeRow(i) = nGroupTypesRow
465	+	endif
466	+	endif
467	+	enddo
468
469	<	if (nGroupTypes.eq.0) then
470	<	nGroupTypes = nGroupTypes + 1
471	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
472	<	groupToGtype(i) = nGroupTypes
469	>	#ifdef IS_MPI
470	>	do j = jstart, jend
471	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
472	>	groupMaxCutoffCol(j) = 0.0_dp
473	>	do ja = groupStartCol(j), groupStartCol(j+1)-1
474	>	atom1 = groupListCol(ja)
475	>
476	>	me_j = atid_col(atom1)
477	>
478	>	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
479	>	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
480	>	endif
481	>	enddo
482	>
483	>	if (nGroupTypesCol.eq.0) then
484	>	nGroupTypesCol = nGroupTypesCol + 1
485	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
486	>	groupToGtypeCol(j) = nGroupTypesCol
487		else
488		GtypeFound = .false.
489	<	do g = 1, nGroupTypes
490	<	if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
491	<	groupToGtype(i) = g
489	>	do g = 1, nGroupTypesCol
490	>	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
491	>	groupToGtypeCol(j) = g
492		GtypeFound = .true.
493		endif
494		enddo
495		if (.not.GtypeFound) then
496	<	nGroupTypes = nGroupTypes + 1
497	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
498	<	groupToGtype(i) = nGroupTypes
496	>	nGroupTypesCol = nGroupTypesCol + 1
497	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
498	>	groupToGtypeCol(j) = nGroupTypesCol
499		endif
500		endif
501		enddo
502
503	+	#else
504	+	! Set pointers to information we just found
505	+	nGroupTypesCol = nGroupTypesRow
506	+	groupToGtypeCol => groupToGtypeRow
507	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
508	+	groupMaxCutoffCol => groupMaxCutoffRow
509	+	#endif
510	+
511		!! allocate the gtypeCutoffMap here.
512	<	allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
512	>	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
513		!! then we do a double loop over all the group TYPES to find the cutoff
514		!! map between groups of two types
515	<
516	<	do i = 1, nGroupTypes
517	<	do j = 1, nGroupTypes
515	>	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
516	>
517	>	do i = 1, nGroupTypesRow
518	>	do j = 1, nGroupTypesCol
519
520		select case(cutoffPolicy)
521		case(TRADITIONAL_CUTOFF_POLICY)
522	<	thisRcut = maxval(gtypeMaxCutoff)
522	>	thisRcut = tradRcut
523		case(MIX_CUTOFF_POLICY)
524	<	thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
524	>	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
525		case(MAX_CUTOFF_POLICY)
526	<	thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
526	>	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
527		case default
528		call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
529		return
530		end select
531		gtypeCutoffMap(i,j)%rcut = thisRcut
532	+
533	+	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
534	+
535		gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
405	–	skin = defaultRlist - defaultRcut
406	–	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
536
537	+	if (.not.haveSkinThickness) then
538	+	skinThickness = 1.0_dp
539	+	endif
540	+
541	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
542	+
543	+	! sanity check
544	+
545	+	if (haveDefaultCutoffs) then
546	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
547	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
548	+	endif
549	+	endif
550		enddo
551		enddo
552	+
553	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
554	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
555	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
556	+	#ifdef IS_MPI
557	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
558	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
559	+	#endif
560	+	groupMaxCutoffCol => null()
561	+	gtypeMaxCutoffCol => null()
562
563		haveGtypeCutoffMap = .true.
564		end subroutine createGtypeCutoffMap
565
566	<	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
415	<	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
416	<	integer, intent(in) :: cutPolicy
566	>	subroutine setCutoffs(defRcut, defRsw, defSP, defSF)
567
568	+	real(kind=dp),intent(in) :: defRcut, defRsw
569	+	logical, intent(in) :: defSP, defSF
570	+	character(len = statusMsgSize) :: errMsg
571	+	integer :: localError
572	+
573		defaultRcut = defRcut
574		defaultRsw = defRsw
575	<	defaultRlist = defRlist
576	<	cutoffPolicy = cutPolicy
577	<	end subroutine setDefaultCutoffs
575	>
576	>	defaultDoShiftPot = defSP
577	>	defaultDoShiftFrc = defSF
578
579	<	subroutine setCutoffPolicy(cutPolicy)
579	>	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
580	>	if (defaultDoShiftFrc) then
581	>	write(errMsg, *) &
582	>	'cutoffRadius and switchingRadius are set to the', newline &
583	>	// tab, 'same value.  OOPSE will use shifted force', newline &
584	>	// tab, 'potentials instead of switching functions.'
585	>
586	>	call handleInfo("setCutoffs", errMsg)
587	>	else
588	>	write(errMsg, *) &
589	>	'cutoffRadius and switchingRadius are set to the', newline &
590	>	// tab, 'same value.  OOPSE will use shifted', newline &
591	>	// tab, 'potentials instead of switching functions.'
592	>
593	>	call handleInfo("setCutoffs", errMsg)
594	>
595	>	defaultDoShiftPot = .true.
596	>	endif
597	>
598	>	endif
599	>
600	>	localError = 0
601	>	call setLJDefaultCutoff( defaultRcut, defaultDoShiftPot, &
602	>	defaultDoShiftFrc )
603	>	call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
604	>	call setCutoffEAM( defaultRcut )
605	>	call setCutoffSC( defaultRcut )
606	>	call set_switch(defaultRsw, defaultRcut)
607	>	call setHmatDangerousRcutValue(defaultRcut)
608	>
609	>	haveDefaultCutoffs = .true.
610	>	haveGtypeCutoffMap = .false.
611
612	+	end subroutine setCutoffs
613	+
614	+	subroutine cWasLame()
615	+
616	+	VisitCutoffsAfterComputing = .true.
617	+	return
618	+
619	+	end subroutine cWasLame
620	+
621	+	subroutine setCutoffPolicy(cutPolicy)
622	+
623		integer, intent(in) :: cutPolicy
624	+
625		cutoffPolicy = cutPolicy
626	<	call createGtypeCutoffMap()
626	>	haveCutoffPolicy = .true.
627	>	haveGtypeCutoffMap = .false.
628	>
629		end subroutine setCutoffPolicy
430	–
630
631	<	subroutine setSimVariables()
433	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
434	<	SIM_uses_EAM = SimUsesEAM()
435	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
436	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
437	<	SIM_uses_PBC = SimUsesPBC()
631	>	subroutine setBoxDipole()
632
633	<	haveSIMvariables = .true.
633	>	do_box_dipole = .true.
634	>
635	>	end subroutine setBoxDipole
636
637	<	return
442	<	end subroutine setSimVariables
637	>	subroutine getBoxDipole( box_dipole )
638
639	+	real(kind=dp), intent(inout), dimension(3) :: box_dipole
640	+
641	+	box_dipole = boxDipole
642	+
643	+	end subroutine getBoxDipole
644	+
645	+	subroutine setElectrostaticMethod( thisESM )
646	+
647	+	integer, intent(in) :: thisESM
648	+
649	+	electrostaticSummationMethod = thisESM
650	+	haveElectrostaticSummationMethod = .true.
651	+
652	+	end subroutine setElectrostaticMethod
653	+
654	+	subroutine setSkinThickness( thisSkin )
655	+
656	+	real(kind=dp), intent(in) :: thisSkin
657	+
658	+	skinThickness = thisSkin
659	+	haveSkinThickness = .true.
660	+	haveGtypeCutoffMap = .false.
661	+
662	+	end subroutine setSkinThickness
663	+
664	+	subroutine setSimVariables()
665	+	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
666	+	SIM_uses_EAM = SimUsesEAM()
667	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
668	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
669	+	SIM_uses_PBC = SimUsesPBC()
670	+	SIM_uses_SC = SimUsesSC()
671	+	SIM_uses_AtomicVirial = SimUsesAtomicVirial()
672	+
673	+	haveSIMvariables = .true.
674	+
675	+	return
676	+	end subroutine setSimVariables
677	+
678		subroutine doReadyCheck(error)
679		integer, intent(out) :: error
446	–
680		integer :: myStatus
681
682		error = 0
683
684		if (.not. haveInteractionHash) then
685	<	myStatus = 0
453	<	call createInteractionHash(myStatus)
454	<	if (myStatus .ne. 0) then
455	<	write(default_error, *) 'createInteractionHash failed in doForces!'
456	<	error = -1
457	<	return
458	<	endif
685	>	call createInteractionHash()
686		endif
687
688		if (.not. haveGtypeCutoffMap) then
689	<	myStatus = 0
690	<	call createGtypeCutoffMap(myStatus)
691	<	if (myStatus .ne. 0) then
692	<	write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
693	<	error = -1
694	<	return
695	<	endif
689	>	call createGtypeCutoffMap()
690	>	endif
691	>
692	>	if (VisitCutoffsAfterComputing) then
693	>	call set_switch(largestRcut, largestRcut)
694	>	call setHmatDangerousRcutValue(largestRcut)
695	>	call setCutoffEAM(largestRcut)
696	>	call setCutoffSC(largestRcut)
697	>	VisitCutoffsAfterComputing = .false.
698		endif
699
700		if (.not. haveSIMvariables) then
701		call setSimVariables()
702		endif
703
475	–	!  if (.not. haveRlist) then
476	–	!     write(default_error, *) 'rList has not been set in doForces!'
477	–	!     error = -1
478	–	!     return
479	–	!  endif
480	–
704		if (.not. haveNeighborList) then
705		write(default_error, *) 'neighbor list has not been initialized in doForces!'
706		error = -1
707		return
708		end if
709	<
709	>
710		if (.not. haveSaneForceField) then
711		write(default_error, *) 'Force Field is not sane in doForces!'
712		error = -1
713		return
714		end if
715	<
715	>
716		#ifdef IS_MPI
717		if (.not. isMPISimSet()) then
718		write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
#	Line 501 | Line 724 | contains
724		end subroutine doReadyCheck
725
726
727	<	subroutine init_FF(thisESM, thisStat)
727	>	subroutine init_FF(thisStat)
728
506	–	integer, intent(in) :: thisESM
507	–	real(kind=dp), intent(in) :: dampingAlpha
729		integer, intent(out) :: thisStat
730		integer :: my_status, nMatches
731		integer, pointer :: MatchList(:) => null()
511	–	real(kind=dp) :: rcut, rrf, rt, dielect
732
733		!! assume things are copacetic, unless they aren't
734		thisStat = 0
735
516	–	electrostaticSummationMethod = thisESM
517	–
736		!! init_FF is called *after* all of the atom types have been
737		!! defined in atype_module using the new_atype subroutine.
738		!!
#	Line 525 | Line 743 | contains
743		FF_uses_Dipoles = .false.
744		FF_uses_GayBerne = .false.
745		FF_uses_EAM = .false.
746	+	FF_uses_SC = .false.
747
748		call getMatchingElementList(atypes, "is_Directional", .true., &
749		nMatches, MatchList)
#	Line 541 | Line 760 | contains
760		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
761		if (nMatches .gt. 0) FF_uses_EAM = .true.
762
763	+	call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
764	+	if (nMatches .gt. 0) FF_uses_SC = .true.
765
766	+
767		haveSaneForceField = .true.
768
547	–	!! check to make sure the reaction field setting makes sense
548	–
549	–	if (FF_uses_Dipoles) then
550	–	if (electrostaticSummationMethod == 3) then
551	–	dielect = getDielect()
552	–	call initialize_rf(dielect)
553	–	endif
554	–	else
555	–	if (electrostaticSummationMethod == 3) then
556	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
557	–	thisStat = -1
558	–	haveSaneForceField = .false.
559	–	return
560	–	endif
561	–	endif
562	–
769		if (FF_uses_EAM) then
770		call init_EAM_FF(my_status)
771		if (my_status /= 0) then
#	Line 570 | Line 776 | contains
776		end if
777		endif
778
573	–	if (FF_uses_GayBerne) then
574	–	call check_gb_pair_FF(my_status)
575	–	if (my_status .ne. 0) then
576	–	thisStat = -1
577	–	haveSaneForceField = .false.
578	–	return
579	–	endif
580	–	endif
581	–
779		if (.not. haveNeighborList) then
780		!! Create neighbor lists
781		call expandNeighborList(nLocal, my_status)
#	Line 612 | Line 809 | contains
809
810		!! Stress Tensor
811		real( kind = dp), dimension(9) :: tau
812	<	real ( kind = dp ) :: pot
812	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
813		logical ( kind = 2) :: do_pot_c, do_stress_c
814		logical :: do_pot
815		logical :: do_stress
816		logical :: in_switching_region
817		#ifdef IS_MPI
818	<	real( kind = DP ) :: pot_local
818	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
819		integer :: nAtomsInRow
820		integer :: nAtomsInCol
821		integer :: nprocs
#	Line 631 | Line 828 | contains
828		integer :: istart, iend
829		integer :: ia, jb, atom1, atom2
830		integer :: nlist
831	<	real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
831	>	real( kind = DP ) :: ratmsq, rgrpsq, rgrp, rag, vpair, vij
832		real( kind = DP ) :: sw, dswdr, swderiv, mf
833	<	real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
834	<	real(kind=dp) :: rfpot, mu_i, virial
833	>	real( kind = DP ) :: rVal
834	>	real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij, fg, dag
835	>	real(kind=dp) :: rfpot, mu_i
836	>	real(kind=dp):: rCut
837		integer :: me_i, me_j, n_in_i, n_in_j
838		logical :: is_dp_i
839		integer :: neighborListSize
#	Line 643 | Line 842 | contains
842		integer :: propPack_i, propPack_j
843		integer :: loopStart, loopEnd, loop
844		integer :: iHash
845	<	real(kind=dp) :: listSkin = 1.0
845	>	integer :: i1
846
847	+	!! the variables for the box dipole moment
848	+	#ifdef IS_MPI
849	+	integer :: pChgCount_local
850	+	integer :: nChgCount_local
851	+	real(kind=dp) :: pChg_local
852	+	real(kind=dp) :: nChg_local
853	+	real(kind=dp), dimension(3) :: pChgPos_local
854	+	real(kind=dp), dimension(3) :: nChgPos_local
855	+	real(kind=dp), dimension(3) :: dipVec_local
856	+	#endif
857	+	integer :: pChgCount
858	+	integer :: nChgCount
859	+	real(kind=dp) :: pChg
860	+	real(kind=dp) :: nChg
861	+	real(kind=dp) :: chg_value
862	+	real(kind=dp), dimension(3) :: pChgPos
863	+	real(kind=dp), dimension(3) :: nChgPos
864	+	real(kind=dp), dimension(3) :: dipVec
865	+	real(kind=dp), dimension(3) :: chgVec
866	+
867	+	!! initialize box dipole variables
868	+	if (do_box_dipole) then
869	+	#ifdef IS_MPI
870	+	pChg_local = 0.0_dp
871	+	nChg_local = 0.0_dp
872	+	pChgCount_local = 0
873	+	nChgCount_local = 0
874	+	do i=1, 3
875	+	pChgPos_local = 0.0_dp
876	+	nChgPos_local = 0.0_dp
877	+	dipVec_local = 0.0_dp
878	+	enddo
879	+	#endif
880	+	pChg = 0.0_dp
881	+	nChg = 0.0_dp
882	+	pChgCount = 0
883	+	nChgCount = 0
884	+	chg_value = 0.0_dp
885	+
886	+	do i=1, 3
887	+	pChgPos(i) = 0.0_dp
888	+	nChgPos(i) = 0.0_dp
889	+	dipVec(i) = 0.0_dp
890	+	chgVec(i) = 0.0_dp
891	+	boxDipole(i) = 0.0_dp
892	+	enddo
893	+	endif
894	+
895		!! initialize local variables
896
897		#ifdef IS_MPI
#	Line 707 | Line 954 | contains
954		! (but only on the first time through):
955		if (loop .eq. loopStart) then
956		#ifdef IS_MPI
957	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
957	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
958		update_nlist)
959		#else
960	<	call checkNeighborList(nGroups, q_group, listSkin, &
960	>	call checkNeighborList(nGroups, q_group, skinThickness, &
961		update_nlist)
962		#endif
963		endif
#	Line 768 | Line 1015 | contains
1015		me_j = atid(j)
1016		call get_interatomic_vector(q_group(:,i), &
1017		q_group(:,j), d_grp, rgrpsq)
1018	<	#endif
1018	>	#endif
1019
1020	<	if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
1020	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
1021		if (update_nlist) then
1022		nlist = nlist + 1
1023
#	Line 790 | Line 1037 | contains
1037
1038		list(nlist) = j
1039		endif
1040	+
1041	+	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
1042
1043	<	if (loop .eq. PAIR_LOOP) then
1044	<	vij = 0.0d0
1045	<	fij(1:3) = 0.0d0
1046	<	endif
1047	<
1048	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
1049	<	in_switching_region)
1050	<
1051	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
1052	<
1053	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
1054	<
1055	<	atom1 = groupListRow(ia)
1056	<
1057	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1058	<
1059	<	atom2 = groupListCol(jb)
1060	<
1061	<	if (skipThisPair(atom1, atom2)) cycle inner
1062	<
1063	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1064	<	d_atm(1:3) = d_grp(1:3)
1065	<	ratmsq = rgrpsq
1066	<	else
1067	<	#ifdef IS_MPI
1068	<	call get_interatomic_vector(q_Row(:,atom1), &
1069	<	q_Col(:,atom2), d_atm, ratmsq)
1043	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
1044	>	if (loop .eq. PAIR_LOOP) then
1045	>	vij = 0.0_dp
1046	>	fij(1) = 0.0_dp
1047	>	fij(2) = 0.0_dp
1048	>	fij(3) = 0.0_dp
1049	>	endif
1050	>
1051	>	call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
1052	>
1053	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
1054	>
1055	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
1056	>
1057	>	atom1 = groupListRow(ia)
1058	>
1059	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1060	>
1061	>	atom2 = groupListCol(jb)
1062	>
1063	>	if (skipThisPair(atom1, atom2))  cycle inner
1064	>
1065	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1066	>	d_atm(1) = d_grp(1)
1067	>	d_atm(2) = d_grp(2)
1068	>	d_atm(3) = d_grp(3)
1069	>	ratmsq = rgrpsq
1070	>	else
1071	>	#ifdef IS_MPI
1072	>	call get_interatomic_vector(q_Row(:,atom1), &
1073	>	q_Col(:,atom2), d_atm, ratmsq)
1074		#else
1075	<	call get_interatomic_vector(q(:,atom1), &
1076	<	q(:,atom2), d_atm, ratmsq)
1075	>	call get_interatomic_vector(q(:,atom1), &
1076	>	q(:,atom2), d_atm, ratmsq)
1077		#endif
1078	<	endif
1079	<
1080	<	if (loop .eq. PREPAIR_LOOP) then
1078	>	endif
1079	>
1080	>	if (loop .eq. PREPAIR_LOOP) then
1081		#ifdef IS_MPI
1082	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1083	<	rgrpsq, d_grp, do_pot, do_stress, &
1084	<	eFrame, A, f, t, pot_local)
1082	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1083	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1084	>	eFrame, A, f, t, pot_local)
1085		#else
1086	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1087	<	rgrpsq, d_grp, do_pot, do_stress, &
1088	<	eFrame, A, f, t, pot)
1086	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1087	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1088	>	eFrame, A, f, t, pot)
1089		#endif
1090	<	else
1090	>	else
1091		#ifdef IS_MPI
1092	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1093	<	do_pot, &
1094	<	eFrame, A, f, t, pot_local, vpair, fpair)
1092	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1093	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1094	>	fpair, d_grp, rgrp, rCut)
1095		#else
1096	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1097	<	do_pot,  &
1098	<	eFrame, A, f, t, pot, vpair, fpair)
1096	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1097	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1098	>	d_grp, rgrp, rCut)
1099		#endif
1100	+	vij = vij + vpair
1101	+	fij(1) = fij(1) + fpair(1)
1102	+	fij(2) = fij(2) + fpair(2)
1103	+	fij(3) = fij(3) + fpair(3)
1104	+	if (do_stress) then
1105	+	call add_stress_tensor(d_atm, fpair, tau)
1106	+	endif
1107	+	endif
1108	+	enddo inner
1109	+	enddo
1110
1111	<	vij = vij + vpair
1112	<	fij(1:3) = fij(1:3) + fpair(1:3)
1113	<	endif
1114	<	enddo inner
1115	<	enddo
1116	<
1117	<	if (loop .eq. PAIR_LOOP) then
1118	<	if (in_switching_region) then
1119	<	swderiv = vij*dswdr/rgrp
1120	<	fij(1) = fij(1) + swderiv*d_grp(1)
1121	<	fij(2) = fij(2) + swderiv*d_grp(2)
1122	<	fij(3) = fij(3) + swderiv*d_grp(3)
1123	<
1124	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
1125	<	atom1=groupListRow(ia)
1126	<	mf = mfactRow(atom1)
1111	>	if (loop .eq. PAIR_LOOP) then
1112	>	if (in_switching_region) then
1113	>	swderiv = vij*dswdr/rgrp
1114	>	fg = swderiv*d_grp
1115	>
1116	>	fij(1) = fij(1) + fg(1)
1117	>	fij(2) = fij(2) + fg(2)
1118	>	fij(3) = fij(3) + fg(3)
1119	>
1120	>	if (do_stress .and. (n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1121	>	call add_stress_tensor(d_atm, fg, tau)
1122	>	endif
1123	>
1124	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1125	>	atom1=groupListRow(ia)
1126	>	mf = mfactRow(atom1)
1127	>	! fg is the force on atom ia due to cutoff group's
1128	>	! presence in switching region
1129	>	fg = swderiv*d_grp*mf
1130		#ifdef IS_MPI
1131	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1132	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1133	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1131	>	f_Row(1,atom1) = f_Row(1,atom1) + fg(1)
1132	>	f_Row(2,atom1) = f_Row(2,atom1) + fg(2)
1133	>	f_Row(3,atom1) = f_Row(3,atom1) + fg(3)
1134		#else
1135	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1136	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1137	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1135	>	f(1,atom1) = f(1,atom1) + fg(1)
1136	>	f(2,atom1) = f(2,atom1) + fg(2)
1137	>	f(3,atom1) = f(3,atom1) + fg(3)
1138		#endif
1139	<	enddo
1140	<
1141	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1142	<	atom2=groupListCol(jb)
877	<	mf = mfactCol(atom2)
1139	>	if (n_in_i .gt. 1) then
1140	>	if (do_stress.and.SIM_uses_AtomicVirial) then
1141	>	! find the distance between the atom and the center of
1142	>	! the cutoff group:
1143		#ifdef IS_MPI
1144	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1145	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
881	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1144	>	call get_interatomic_vector(q_Row(:,atom1), &
1145	>	q_group_Row(:,i), dag, rag)
1146		#else
1147	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1148	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
885	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1147	>	call get_interatomic_vector(q(:,atom1), &
1148	>	q_group(:,i), dag, rag)
1149		#endif
1150	<	enddo
1150	>	call add_stress_tensor(dag,fg,tau)
1151	>	endif
1152	>	endif
1153	>	enddo
1154	>
1155	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1156	>	atom2=groupListCol(jb)
1157	>	mf = mfactCol(atom2)
1158	>	! fg is the force on atom jb due to cutoff group's
1159	>	! presence in switching region
1160	>	fg = -swderiv*d_grp*mf
1161	>	#ifdef IS_MPI
1162	>	f_Col(1,atom2) = f_Col(1,atom2) + fg(1)
1163	>	f_Col(2,atom2) = f_Col(2,atom2) + fg(2)
1164	>	f_Col(3,atom2) = f_Col(3,atom2) + fg(3)
1165	>	#else
1166	>	f(1,atom2) = f(1,atom2) + fg(1)
1167	>	f(2,atom2) = f(2,atom2) + fg(2)
1168	>	f(3,atom2) = f(3,atom2) + fg(3)
1169	>	#endif
1170	>	if (n_in_j .gt. 1) then
1171	>	if (do_stress.and.SIM_uses_AtomicVirial) then
1172	>	! find the distance between the atom and the center of
1173	>	! the cutoff group:
1174	>	#ifdef IS_MPI
1175	>	call get_interatomic_vector(q_Col(:,atom2), &
1176	>	q_group_Col(:,j), dag, rag)
1177	>	#else
1178	>	call get_interatomic_vector(q(:,atom2), &
1179	>	q_group(:,j), dag, rag)
1180	>	#endif
1181	>	call add_stress_tensor(dag,fg,tau)
1182	>	endif
1183	>	endif
1184	>	enddo
1185	>	endif
1186		endif
889	–
890	–	if (do_stress) call add_stress_tensor(d_grp, fij)
1187		endif
1188	<	end if
1188	>	endif
1189		enddo
1190	+
1191		enddo outer
1192
1193		if (update_nlist) then
#	Line 908 | Line 1205 | contains
1205		endif
1206
1207		if (loop .eq. PREPAIR_LOOP) then
1208	+	#ifdef IS_MPI
1209	+	call do_preforce(nlocal, pot_local)
1210	+	#else
1211		call do_preforce(nlocal, pot)
1212	+	#endif
1213		endif
1214
1215		enddo
#	Line 950 | Line 1251 | contains
1251
1252		if (do_pot) then
1253		! scatter/gather pot_row into the members of my column
1254	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1255	<
1254	>	do i = 1,LR_POT_TYPES
1255	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1256	>	end do
1257		! scatter/gather pot_local into all other procs
1258		! add resultant to get total pot
1259		do i = 1, nlocal
1260	<	pot_local = pot_local + pot_Temp(i)
1260	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1261	>	+ pot_Temp(1:LR_POT_TYPES,i)
1262		enddo
1263
1264		pot_Temp = 0.0_DP
1265	<
1266	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1265	>	do i = 1,LR_POT_TYPES
1266	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1267	>	end do
1268		do i = 1, nlocal
1269	<	pot_local = pot_local + pot_Temp(i)
1269	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1270	>	+ pot_Temp(1:LR_POT_TYPES,i)
1271		enddo
1272
1273		endif
1274		#endif
1275
1276	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1276	>	if (SIM_requires_postpair_calc) then
1277	>	do i = 1, nlocal
1278	>
1279	>	! we loop only over the local atoms, so we don't need row and column
1280	>	! lookups for the types
1281	>
1282	>	me_i = atid(i)
1283	>
1284	>	! is the atom electrostatic?  See if it would have an
1285	>	! electrostatic interaction with itself
1286	>	iHash = InteractionHash(me_i,me_i)
1287
1288	<	if (electrostaticSummationMethod == 3) then
974	<
1288	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1289		#ifdef IS_MPI
1290	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1291	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1292	<	do i = 1,nlocal
1293	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1294	<	end do
1290	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1291	>	t, do_pot)
1292	>	#else
1293	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1294	>	t, do_pot)
1295		#endif
1296	<
1297	<	do i = 1, nLocal
1298	<
1299	<	rfpot = 0.0_DP
1296	>	endif
1297	>
1298	>
1299	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1300	>
1301	>	! loop over the excludes to accumulate RF stuff we've
1302	>	! left out of the normal pair loop
1303	>
1304	>	do i1 = 1, nSkipsForAtom(i)
1305	>	j = skipsForAtom(i, i1)
1306	>
1307	>	! prevent overcounting of the skips
1308	>	if (i.lt.j) then
1309	>	call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1310	>	rVal = sqrt(ratmsq)
1311	>	call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1312		#ifdef IS_MPI
1313	<	me_i = atid_row(i)
1313	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1314	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1315		#else
1316	<	me_i = atid(i)
1316	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1317	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1318		#endif
1319	<	iHash = InteractionHash(me_i,me_j)
1320	<
1321	<	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1319	>	endif
1320	>	enddo
1321	>	endif
1322
1323	<	mu_i = getDipoleMoment(me_i)
996	<
997	<	!! The reaction field needs to include a self contribution
998	<	!! to the field:
999	<	call accumulate_self_rf(i, mu_i, eFrame)
1000	<	!! Get the reaction field contribution to the
1001	<	!! potential and torques:
1002	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1323	>	if (do_box_dipole) then
1324		#ifdef IS_MPI
1325	<	pot_local = pot_local + rfpot
1325	>	call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1326	>	nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1327	>	pChgCount_local, nChgCount_local)
1328		#else
1329	<	pot = pot + rfpot
1330	<
1329	>	call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1330	>	pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1331		#endif
1332	<	endif
1333	<	enddo
1011	<	endif
1332	>	endif
1333	>	enddo
1334		endif
1335
1014	–
1336		#ifdef IS_MPI
1016	–
1337		if (do_pot) then
1338	<	pot = pot + pot_local
1339	<	!! we assume the c code will do the allreduce to get the total potential
1340	<	!! we could do it right here if we needed to...
1338	>	#ifdef SINGLE_PRECISION
1339	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1340	>	mpi_comm_world,mpi_err)
1341	>	#else
1342	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1343	>	mpi_sum, mpi_comm_world,mpi_err)
1344	>	#endif
1345		endif
1346	+
1347	+	if (do_box_dipole) then
1348
1349	<	if (do_stress) then
1350	<	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1351	<	mpi_comm_world,mpi_err)
1352	<	call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1353	<	mpi_comm_world,mpi_err)
1354	<	endif
1355	<
1349	>	#ifdef SINGLE_PRECISION
1350	>	call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1351	>	mpi_comm_world, mpi_err)
1352	>	call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1353	>	mpi_comm_world, mpi_err)
1354	>	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1355	>	mpi_comm_world, mpi_err)
1356	>	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1357	>	mpi_comm_world, mpi_err)
1358	>	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1359	>	mpi_comm_world, mpi_err)
1360	>	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1361	>	mpi_comm_world, mpi_err)
1362	>	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1363	>	mpi_comm_world, mpi_err)
1364		#else
1365	+	call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1366	+	mpi_comm_world, mpi_err)
1367	+	call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1368	+	mpi_comm_world, mpi_err)
1369	+	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1370	+	mpi_sum, mpi_comm_world, mpi_err)
1371	+	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1372	+	mpi_sum, mpi_comm_world, mpi_err)
1373	+	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1374	+	mpi_sum, mpi_comm_world, mpi_err)
1375	+	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1376	+	mpi_sum, mpi_comm_world, mpi_err)
1377	+	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1378	+	mpi_sum, mpi_comm_world, mpi_err)
1379	+	#endif
1380
1032	–	if (do_stress) then
1033	–	tau = tau_Temp
1034	–	virial = virial_Temp
1381		endif
1382	<
1382	>
1383		#endif
1384
1385	+	if (do_box_dipole) then
1386	+	! first load the accumulated dipole moment (if dipoles were present)
1387	+	boxDipole(1) = dipVec(1)
1388	+	boxDipole(2) = dipVec(2)
1389	+	boxDipole(3) = dipVec(3)
1390	+
1391	+	! now include the dipole moment due to charges
1392	+	! use the lesser of the positive and negative charge totals
1393	+	if (nChg .le. pChg) then
1394	+	chg_value = nChg
1395	+	else
1396	+	chg_value = pChg
1397	+	endif
1398	+
1399	+	! find the average positions
1400	+	if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1401	+	pChgPos = pChgPos / pChgCount
1402	+	nChgPos = nChgPos / nChgCount
1403	+	endif
1404	+
1405	+	! dipole is from the negative to the positive (physics notation)
1406	+	chgVec(1) = pChgPos(1) - nChgPos(1)
1407	+	chgVec(2) = pChgPos(2) - nChgPos(2)
1408	+	chgVec(3) = pChgPos(3) - nChgPos(3)
1409	+
1410	+	boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1411	+	boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1412	+	boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1413	+
1414	+	endif
1415	+
1416		end subroutine do_force_loop
1417
1418		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1419	<	eFrame, A, f, t, pot, vpair, fpair)
1419	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1420
1421	<	real( kind = dp ) :: pot, vpair, sw
1421	>	real( kind = dp ) :: vpair, sw
1422	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1423		real( kind = dp ), dimension(3) :: fpair
1424		real( kind = dp ), dimension(nLocal)   :: mfact
1425		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1052 | Line 1430 | contains
1430		logical, intent(inout) :: do_pot
1431		integer, intent(in) :: i, j
1432		real ( kind = dp ), intent(inout) :: rijsq
1433	<	real ( kind = dp )                :: r
1433	>	real ( kind = dp ), intent(inout) :: r_grp
1434		real ( kind = dp ), intent(inout) :: d(3)
1435	+	real ( kind = dp ), intent(inout) :: d_grp(3)
1436	+	real ( kind = dp ), intent(inout) :: rCut
1437	+	real ( kind = dp ) :: r
1438	+	real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1439		integer :: me_i, me_j
1440	+	integer :: k
1441
1442		integer :: iHash
1443
1444		r = sqrt(rijsq)
1445	<	vpair = 0.0d0
1446	<	fpair(1:3) = 0.0d0
1445	>
1446	>	vpair = 0.0_dp
1447	>	fpair(1:3) = 0.0_dp
1448
1449		#ifdef IS_MPI
1450		me_i = atid_row(i)
#	Line 1071 | Line 1455 | contains
1455		#endif
1456
1457		iHash = InteractionHash(me_i, me_j)
1458	<
1458	>
1459		if ( iand(iHash, LJ_PAIR).ne.0 ) then
1460	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1460	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1461	>	pot(VDW_POT), f, do_pot)
1462		endif
1463	<
1463	>
1464		if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1465	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1466	<	pot, eFrame, f, t, do_pot)
1082	<
1083	<	if (electrostaticSummationMethod == 3) then
1084	<
1085	<	! CHECK ME (RF needs to know about all electrostatic types)
1086	<	call accumulate_rf(i, j, r, eFrame, sw)
1087	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1088	<	endif
1089	<
1465	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1466	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1467		endif
1468	<
1468	>
1469		if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1470		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1471	<	pot, A, f, t, do_pot)
1471	>	pot(HB_POT), A, f, t, do_pot)
1472		endif
1473	<
1473	>
1474		if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1475		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1476	<	pot, A, f, t, do_pot)
1476	>	pot(HB_POT), A, f, t, do_pot)
1477		endif
1478	<
1478	>
1479		if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1480		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1481	<	pot, A, f, t, do_pot)
1481	>	pot(VDW_POT), A, f, t, do_pot)
1482		endif
1483
1484		if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1485	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1486	<	!           pot, A, f, t, do_pot)
1485	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1486	>	pot(VDW_POT), A, f, t, do_pot)
1487		endif
1488	<
1488	>
1489		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1490	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1491	<	do_pot)
1490	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1491	>	pot(METALLIC_POT), f, do_pot)
1492		endif
1493	<
1493	>
1494		if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1495		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1496	<	pot, A, f, t, do_pot)
1496	>	pot(VDW_POT), A, f, t, do_pot)
1497		endif
1498	<
1498	>
1499		if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1500		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1501	<	pot, A, f, t, do_pot)
1501	>	pot(VDW_POT), A, f, t, do_pot)
1502		endif
1503	<
1503	>
1504	>	if ( iand(iHash, SC_PAIR).ne.0 ) then
1505	>	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1506	>	pot(METALLIC_POT), f, do_pot)
1507	>	endif
1508	>
1509		end subroutine do_pair
1510
1511	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1511	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1512		do_pot, do_stress, eFrame, A, f, t, pot)
1513
1514	<	real( kind = dp ) :: pot, sw
1514	>	real( kind = dp ) :: sw
1515	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1516		real( kind = dp ), dimension(9,nLocal) :: eFrame
1517		real (kind=dp), dimension(9,nLocal) :: A
1518		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1137 | Line 1520 | contains
1520
1521		logical, intent(inout) :: do_pot, do_stress
1522		integer, intent(in) :: i, j
1523	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1523	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1524		real ( kind = dp )                :: r, rc
1525		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1526
1527		integer :: me_i, me_j, iHash
1528
1529		r = sqrt(rijsq)
1530	<
1530	>
1531		#ifdef IS_MPI
1532		me_i = atid_row(i)
1533		me_j = atid_col(j)
#	Line 1156 | Line 1539 | contains
1539		iHash = InteractionHash(me_i, me_j)
1540
1541		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1542	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1542	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1543		endif
1544	+
1545	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1546	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1547	+	endif
1548
1549		end subroutine do_prepair
1550
1551
1552		subroutine do_preforce(nlocal,pot)
1553		integer :: nlocal
1554	<	real( kind = dp ) :: pot
1554	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1555
1556		if (FF_uses_EAM .and. SIM_uses_EAM) then
1557	<	call calc_EAM_preforce_Frho(nlocal,pot)
1557	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1558		endif
1559	<
1560	<
1559	>	if (FF_uses_SC .and. SIM_uses_SC) then
1560	>	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1561	>	endif
1562		end subroutine do_preforce
1563
1564
#	Line 1182 | Line 1570 | contains
1570		real( kind = dp ) :: d(3), scaled(3)
1571		integer i
1572
1573	<	d(1:3) = q_j(1:3) - q_i(1:3)
1573	>	d(1) = q_j(1) - q_i(1)
1574	>	d(2) = q_j(2) - q_i(2)
1575	>	d(3) = q_j(3) - q_i(3)
1576
1577		! Wrap back into periodic box if necessary
1578		if ( SIM_uses_PBC ) then
1579
1580		if( .not.boxIsOrthorhombic ) then
1581		! calc the scaled coordinates.
1582	+	! scaled = matmul(HmatInv, d)
1583
1584	<	scaled = matmul(HmatInv, d)
1585	<
1584	>	scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1585	>	scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1586	>	scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1587	>
1588		! wrap the scaled coordinates
1589
1590	<	scaled = scaled  - anint(scaled)
1590	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1591	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1592	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1593
1199	–
1594		! calc the wrapped real coordinates from the wrapped scaled
1595		! coordinates
1596	+	! d = matmul(Hmat,scaled)
1597	+	d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1598	+	d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1599	+	d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1600
1203	–	d = matmul(Hmat,scaled)
1204	–
1601		else
1602		! calc the scaled coordinates.
1603
1604	<	do i = 1, 3
1605	<	scaled(i) = d(i) * HmatInv(i,i)
1606	<
1607	<	! wrap the scaled coordinates
1604	>	scaled(1) = d(1) * HmatInv(1,1)
1605	>	scaled(2) = d(2) * HmatInv(2,2)
1606	>	scaled(3) = d(3) * HmatInv(3,3)
1607	>
1608	>	! wrap the scaled coordinates
1609	>
1610	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1611	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1612	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1613
1614	<	scaled(i) = scaled(i) - anint(scaled(i))
1614	>	! calc the wrapped real coordinates from the wrapped scaled
1615	>	! coordinates
1616
1617	<	! calc the wrapped real coordinates from the wrapped scaled
1618	<	! coordinates
1617	>	d(1) = scaled(1)*Hmat(1,1)
1618	>	d(2) = scaled(2)*Hmat(2,2)
1619	>	d(3) = scaled(3)*Hmat(3,3)
1620
1218	–	d(i) = scaled(i)*Hmat(i,i)
1219	–	enddo
1621		endif
1622
1623		endif
1624
1625	<	r_sq = dot_product(d,d)
1625	>	r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1626
1627		end subroutine get_interatomic_vector
1628
#	Line 1253 | Line 1654 | contains
1654		pot_Col = 0.0_dp
1655		pot_Temp = 0.0_dp
1656
1256	–	rf_Row = 0.0_dp
1257	–	rf_Col = 0.0_dp
1258	–	rf_Temp = 0.0_dp
1259	–
1657		#endif
1658
1659		if (FF_uses_EAM .and. SIM_uses_EAM) then
1660		call clean_EAM()
1661		endif
1662
1266	–	rf = 0.0_dp
1267	–	tau_Temp = 0.0_dp
1268	–	virial_Temp = 0.0_dp
1663		end subroutine zero_work_arrays
1664
1665		function skipThisPair(atom1, atom2) result(skip_it)
#	Line 1357 | Line 1751 | contains
1751
1752		function FF_RequiresPrepairCalc() result(doesit)
1753		logical :: doesit
1754	<	doesit = FF_uses_EAM
1754	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1755	>	.or. FF_uses_MEAM
1756		end function FF_RequiresPrepairCalc
1757
1363	–	function FF_RequiresPostpairCalc() result(doesit)
1364	–	logical :: doesit
1365	–	if (electrostaticSummationMethod == 3) doesit = .true.
1366	–	end function FF_RequiresPostpairCalc
1367	–
1758		#ifdef PROFILE
1759		function getforcetime() result(totalforcetime)
1760		real(kind=dp) :: totalforcetime
#	Line 1374 | Line 1764 | contains
1764
1765		!! This cleans componets of force arrays belonging only to fortran
1766
1767	<	subroutine add_stress_tensor(dpair, fpair)
1767	>	subroutine add_stress_tensor(dpair, fpair, tau)
1768
1769		real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1770	+	real( kind = dp ), dimension(9), intent(inout) :: tau
1771
1772		! because the d vector is the rj - ri vector, and
1773		! because fx, fy, fz are the force on atom i, we need a
1774		! negative sign here:
1775
1776	<	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1777	<	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1778	<	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1779	<	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1780	<	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1781	<	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1782	<	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1783	<	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1784	<	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1776	>	tau(1) = tau(1) - dpair(1) * fpair(1)
1777	>	tau(2) = tau(2) - dpair(1) * fpair(2)
1778	>	tau(3) = tau(3) - dpair(1) * fpair(3)
1779	>	tau(4) = tau(4) - dpair(2) * fpair(1)
1780	>	tau(5) = tau(5) - dpair(2) * fpair(2)
1781	>	tau(6) = tau(6) - dpair(2) * fpair(3)
1782	>	tau(7) = tau(7) - dpair(3) * fpair(1)
1783	>	tau(8) = tau(8) - dpair(3) * fpair(2)
1784	>	tau(9) = tau(9) - dpair(3) * fpair(3)
1785
1395	–	virial_Temp = virial_Temp + &
1396	–	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1397	–
1786		end subroutine add_stress_tensor
1787
1788		end module doForces

Diff Legend

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