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 3127 by gezelter, Mon Apr 9 18:24:00 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.86 2007-04-09 18:24:00 gezelter Exp $, $Date: 2007-04-09 18:24:00 $, $Name: not supported by cvs2svn $, $Revision: 1.86 $
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 :: defaultDoShift
117	+
118		public :: init_FF
119	<	public :: setDefaultCutoffs
119	>	public :: setCutoffs
120	>	public :: cWasLame
121	>	public :: setElectrostaticMethod
122	>	public :: setBoxDipole
123	>	public :: getBoxDipole
124	>	public :: setCutoffPolicy
125	>	public :: setSkinThickness
126		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
127
128		#ifdef PROFILE
129		public :: getforcetime
#	Line 122 | Line 136 | module doForces
136		! Bit hash to determine pair-pair interactions.
137		integer, dimension(:,:), allocatable :: InteractionHash
138		real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
139	<	real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
140	<	integer, dimension(:), allocatable :: groupToGtype
141	<	real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
139	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
140	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
141	>
142	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
143	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
144	>
145	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
146	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
147		type ::gtypeCutoffs
148		real(kind=dp) :: rcut
149		real(kind=dp) :: rcutsq
#	Line 132 | Line 151 | module doForces
151		end type gtypeCutoffs
152		type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
153
154	<	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
155	<	real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
137	<
154	>	real(kind=dp), dimension(3) :: boxDipole
155	>
156		contains
157
158	<	subroutine createInteractionHash(status)
158	>	subroutine createInteractionHash()
159		integer :: nAtypes
142	–	integer, intent(out) :: status
160		integer :: i
161		integer :: j
162		integer :: iHash
#	Line 151 | Line 168 | contains
168		logical :: i_is_GB
169		logical :: i_is_EAM
170		logical :: i_is_Shape
171	+	logical :: i_is_SC
172	+	logical :: i_is_MEAM
173		logical :: j_is_LJ
174		logical :: j_is_Elect
175		logical :: j_is_Sticky
#	Line 158 | Line 177 | contains
177		logical :: j_is_GB
178		logical :: j_is_EAM
179		logical :: j_is_Shape
180	+	logical :: j_is_SC
181	+	logical :: j_is_MEAM
182		real(kind=dp) :: myRcut
183
163	–	status = 0
164	–
184		if (.not. associated(atypes)) then
185	<	call handleError("atype", "atypes was not present before call of createInteractionHash!")
167	<	status = -1
185	>	call handleError("doForces", "atypes was not present before call of createInteractionHash!")
186		return
187		endif
188
189		nAtypes = getSize(atypes)
190
191		if (nAtypes == 0) then
192	<	status = -1
192	>	call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
193		return
194		end if
195
196		if (.not. allocated(InteractionHash)) then
197		allocate(InteractionHash(nAtypes,nAtypes))
198	+	else
199	+	deallocate(InteractionHash)
200	+	allocate(InteractionHash(nAtypes,nAtypes))
201		endif
202
203		if (.not. allocated(atypeMaxCutoff)) then
204		allocate(atypeMaxCutoff(nAtypes))
205	+	else
206	+	deallocate(atypeMaxCutoff)
207	+	allocate(atypeMaxCutoff(nAtypes))
208		endif
209
210		do i = 1, nAtypes
#	Line 191 | Line 215 | contains
215		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
216		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
217		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
218	+	call getElementProperty(atypes, i, "is_SC", i_is_SC)
219	+	call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
220
221		do j = i, nAtypes
222
#	Line 204 | Line 230 | contains
230		call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
231		call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
232		call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
233	+	call getElementProperty(atypes, j, "is_SC", j_is_SC)
234	+	call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
235
236		if (i_is_LJ .and. j_is_LJ) then
237		iHash = ior(iHash, LJ_PAIR)
#	Line 223 | Line 251 | contains
251
252		if (i_is_EAM .and. j_is_EAM) then
253		iHash = ior(iHash, EAM_PAIR)
254	+	endif
255	+
256	+	if (i_is_SC .and. j_is_SC) then
257	+	iHash = ior(iHash, SC_PAIR)
258		endif
259
260		if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
#	Line 244 | Line 276 | contains
276		haveInteractionHash = .true.
277		end subroutine createInteractionHash
278
279	<	subroutine createGtypeCutoffMap(stat)
279	>	subroutine createGtypeCutoffMap()
280
249	–	integer, intent(out), optional :: stat
281		logical :: i_is_LJ
282		logical :: i_is_Elect
283		logical :: i_is_Sticky
#	Line 254 | Line 285 | contains
285		logical :: i_is_GB
286		logical :: i_is_EAM
287		logical :: i_is_Shape
288	+	logical :: i_is_SC
289		logical :: GtypeFound
290
291		integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
292	<	integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
292	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
293		integer :: nGroupsInRow
294	<	real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
294	>	integer :: nGroupsInCol
295	>	integer :: nGroupTypesRow,nGroupTypesCol
296	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
297		real(kind=dp) :: biggestAtypeCutoff
298
265	–	stat = 0
299		if (.not. haveInteractionHash) then
300	<	call createInteractionHash(myStatus)
268	<	if (myStatus .ne. 0) then
269	<	write(default_error, *) 'createInteractionHash failed in doForces!'
270	<	stat = -1
271	<	return
272	<	endif
300	>	call createInteractionHash()
301		endif
302		#ifdef IS_MPI
303		nGroupsInRow = getNgroupsInRow(plan_group_row)
304	+	nGroupsInCol = getNgroupsInCol(plan_group_col)
305		#endif
306		nAtypes = getSize(atypes)
307		! Set all of the initial cutoffs to zero.
#	Line 286 | Line 315 | contains
315		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
316		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
317		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
318	<
318	>	call getElementProperty(atypes, i, "is_SC", i_is_SC)
319
320	<	if (i_is_LJ) then
321	<	thisRcut = getSigma(i) * 2.5_dp
322	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
323	<	endif
324	<	if (i_is_Elect) then
325	<	thisRcut = defaultRcut
326	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
327	<	endif
328	<	if (i_is_Sticky) then
329	<	thisRcut = getStickyCut(i)
330	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
331	<	endif
332	<	if (i_is_StickyP) then
333	<	thisRcut = getStickyPowerCut(i)
334	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
335	<	endif
336	<	if (i_is_GB) then
337	<	thisRcut = getGayBerneCut(i)
338	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
339	<	endif
340	<	if (i_is_EAM) then
341	<	thisRcut = getEAMCut(i)
342	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
343	<	endif
344	<	if (i_is_Shape) then
345	<	thisRcut = getShapeCut(i)
346	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
320	>	if (haveDefaultCutoffs) then
321	>	atypeMaxCutoff(i) = defaultRcut
322	>	else
323	>	if (i_is_LJ) then
324	>	thisRcut = getSigma(i) * 2.5_dp
325	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326	>	endif
327	>	if (i_is_Elect) then
328	>	thisRcut = defaultRcut
329	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330	>	endif
331	>	if (i_is_Sticky) then
332	>	thisRcut = getStickyCut(i)
333	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334	>	endif
335	>	if (i_is_StickyP) then
336	>	thisRcut = getStickyPowerCut(i)
337	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
338	>	endif
339	>	if (i_is_GB) then
340	>	thisRcut = getGayBerneCut(i)
341	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
342	>	endif
343	>	if (i_is_EAM) then
344	>	thisRcut = getEAMCut(i)
345	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
346	>	endif
347	>	if (i_is_Shape) then
348	>	thisRcut = getShapeCut(i)
349	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
350	>	endif
351	>	if (i_is_SC) then
352	>	thisRcut = getSCCut(i)
353	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
354	>	endif
355		endif
356	<
356	>
357		if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
358		biggestAtypeCutoff = atypeMaxCutoff(i)
359		endif
360	+
361		endif
362		enddo
325	–
326	–	nGroupTypes = 0
363
364		istart = 1
365	+	jstart = 1
366		#ifdef IS_MPI
367		iend = nGroupsInRow
368	+	jend = nGroupsInCol
369		#else
370		iend = nGroups
371	+	jend = nGroups
372		#endif
373
374		!! allocate the groupToGtype and gtypeMaxCutoff here.
375	<	if(.not.allocated(groupToGtype)) then
376	<	allocate(groupToGtype(iend))
377	<	allocate(groupMaxCutoff(iend))
378	<	allocate(gtypeMaxCutoff(iend))
379	<	groupMaxCutoff = 0.0_dp
380	<	gtypeMaxCutoff = 0.0_dp
375	>	if(.not.allocated(groupToGtypeRow)) then
376	>	!  allocate(groupToGtype(iend))
377	>	allocate(groupToGtypeRow(iend))
378	>	else
379	>	deallocate(groupToGtypeRow)
380	>	allocate(groupToGtypeRow(iend))
381		endif
382	+	if(.not.allocated(groupMaxCutoffRow)) then
383	+	allocate(groupMaxCutoffRow(iend))
384	+	else
385	+	deallocate(groupMaxCutoffRow)
386	+	allocate(groupMaxCutoffRow(iend))
387	+	end if
388	+
389	+	if(.not.allocated(gtypeMaxCutoffRow)) then
390	+	allocate(gtypeMaxCutoffRow(iend))
391	+	else
392	+	deallocate(gtypeMaxCutoffRow)
393	+	allocate(gtypeMaxCutoffRow(iend))
394	+	endif
395	+
396	+
397	+	#ifdef IS_MPI
398	+	! We only allocate new storage if we are in MPI because Ncol /= Nrow
399	+	if(.not.associated(groupToGtypeCol)) then
400	+	allocate(groupToGtypeCol(jend))
401	+	else
402	+	deallocate(groupToGtypeCol)
403	+	allocate(groupToGtypeCol(jend))
404	+	end if
405	+
406	+	if(.not.associated(groupMaxCutoffCol)) then
407	+	allocate(groupMaxCutoffCol(jend))
408	+	else
409	+	deallocate(groupMaxCutoffCol)
410	+	allocate(groupMaxCutoffCol(jend))
411	+	end if
412	+	if(.not.associated(gtypeMaxCutoffCol)) then
413	+	allocate(gtypeMaxCutoffCol(jend))
414	+	else
415	+	deallocate(gtypeMaxCutoffCol)
416	+	allocate(gtypeMaxCutoffCol(jend))
417	+	end if
418	+
419	+	groupMaxCutoffCol = 0.0_dp
420	+	gtypeMaxCutoffCol = 0.0_dp
421	+
422	+	#endif
423	+	groupMaxCutoffRow = 0.0_dp
424	+	gtypeMaxCutoffRow = 0.0_dp
425	+
426	+
427		!! first we do a single loop over the cutoff groups to find the
428		!! largest cutoff for any atypes present in this group.  We also
429		!! create gtypes at this point.
430
431	<	tol = 1.0d-6
432	<
431	>	tol = 1.0e-6_dp
432	>	nGroupTypesRow = 0
433	>	nGroupTypesCol = 0
434		do i = istart, iend
435		n_in_i = groupStartRow(i+1) - groupStartRow(i)
436	<	groupMaxCutoff(i) = 0.0_dp
436	>	groupMaxCutoffRow(i) = 0.0_dp
437		do ia = groupStartRow(i), groupStartRow(i+1)-1
438		atom1 = groupListRow(ia)
439		#ifdef IS_MPI
#	Line 356 | Line 441 | contains
441		#else
442		me_i = atid(atom1)
443		#endif
444	<	if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
445	<	groupMaxCutoff(i)=atypeMaxCutoff(me_i)
444	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
445	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
446		endif
447		enddo
448	+	if (nGroupTypesRow.eq.0) then
449	+	nGroupTypesRow = nGroupTypesRow + 1
450	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
451	+	groupToGtypeRow(i) = nGroupTypesRow
452	+	else
453	+	GtypeFound = .false.
454	+	do g = 1, nGroupTypesRow
455	+	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
456	+	groupToGtypeRow(i) = g
457	+	GtypeFound = .true.
458	+	endif
459	+	enddo
460	+	if (.not.GtypeFound) then
461	+	nGroupTypesRow = nGroupTypesRow + 1
462	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
463	+	groupToGtypeRow(i) = nGroupTypesRow
464	+	endif
465	+	endif
466	+	enddo
467
468	<	if (nGroupTypes.eq.0) then
469	<	nGroupTypes = nGroupTypes + 1
470	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
471	<	groupToGtype(i) = nGroupTypes
468	>	#ifdef IS_MPI
469	>	do j = jstart, jend
470	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
471	>	groupMaxCutoffCol(j) = 0.0_dp
472	>	do ja = groupStartCol(j), groupStartCol(j+1)-1
473	>	atom1 = groupListCol(ja)
474	>
475	>	me_j = atid_col(atom1)
476	>
477	>	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
478	>	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
479	>	endif
480	>	enddo
481	>
482	>	if (nGroupTypesCol.eq.0) then
483	>	nGroupTypesCol = nGroupTypesCol + 1
484	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
485	>	groupToGtypeCol(j) = nGroupTypesCol
486		else
487		GtypeFound = .false.
488	<	do g = 1, nGroupTypes
489	<	if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
490	<	groupToGtype(i) = g
488	>	do g = 1, nGroupTypesCol
489	>	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
490	>	groupToGtypeCol(j) = g
491		GtypeFound = .true.
492		endif
493		enddo
494		if (.not.GtypeFound) then
495	<	nGroupTypes = nGroupTypes + 1
496	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
497	<	groupToGtype(i) = nGroupTypes
495	>	nGroupTypesCol = nGroupTypesCol + 1
496	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
497	>	groupToGtypeCol(j) = nGroupTypesCol
498		endif
499		endif
500		enddo
501
502	+	#else
503	+	! Set pointers to information we just found
504	+	nGroupTypesCol = nGroupTypesRow
505	+	groupToGtypeCol => groupToGtypeRow
506	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
507	+	groupMaxCutoffCol => groupMaxCutoffRow
508	+	#endif
509	+
510		!! allocate the gtypeCutoffMap here.
511	<	allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
511	>	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
512		!! then we do a double loop over all the group TYPES to find the cutoff
513		!! map between groups of two types
514	<
515	<	do i = 1, nGroupTypes
516	<	do j = 1, nGroupTypes
514	>	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
515	>
516	>	do i = 1, nGroupTypesRow
517	>	do j = 1, nGroupTypesCol
518
519		select case(cutoffPolicy)
520		case(TRADITIONAL_CUTOFF_POLICY)
521	<	thisRcut = maxval(gtypeMaxCutoff)
521	>	thisRcut = tradRcut
522		case(MIX_CUTOFF_POLICY)
523	<	thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
523	>	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
524		case(MAX_CUTOFF_POLICY)
525	<	thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
525	>	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
526		case default
527		call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
528		return
529		end select
530		gtypeCutoffMap(i,j)%rcut = thisRcut
531	<	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
532	<	skin = defaultRlist - defaultRcut
406	<	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
531	>
532	>	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
533
534	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
535	+
536	+	if (.not.haveSkinThickness) then
537	+	skinThickness = 1.0_dp
538	+	endif
539	+
540	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
541	+
542	+	! sanity check
543	+
544	+	if (haveDefaultCutoffs) then
545	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
546	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
547	+	endif
548	+	endif
549		enddo
550		enddo
551	+
552	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
553	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
554	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
555	+	#ifdef IS_MPI
556	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
557	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
558	+	#endif
559	+	groupMaxCutoffCol => null()
560	+	gtypeMaxCutoffCol => null()
561
562		haveGtypeCutoffMap = .true.
563		end subroutine createGtypeCutoffMap
564
565	<	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
415	<	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
416	<	integer, intent(in) :: cutPolicy
565	>	subroutine setCutoffs(defRcut, defRsw)
566
567	+	real(kind=dp),intent(in) :: defRcut, defRsw
568	+	character(len = statusMsgSize) :: errMsg
569	+	integer :: localError
570	+
571		defaultRcut = defRcut
572		defaultRsw = defRsw
573	<	defaultRlist = defRlist
574	<	cutoffPolicy = cutPolicy
575	<	end subroutine setDefaultCutoffs
573	>
574	>	defaultDoShift = .false.
575	>	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
576	>
577	>	write(errMsg, *) &
578	>	'cutoffRadius and switchingRadius are set to the same', newline &
579	>	// tab, 'value.  OOPSE will use shifted force van der Waals', newline &
580	>	// tab, 'potentials instead of switching functions.'
581	>
582	>	call handleInfo("setCutoffs", errMsg)
583	>
584	>	defaultDoShift = .true.
585	>
586	>	endif
587	>
588	>	localError = 0
589	>	call setLJDefaultCutoff( defaultRcut, defaultDoShift )
590	>	call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
591	>	call setCutoffEAM( defaultRcut )
592	>	call setCutoffSC( defaultRcut )
593	>	call set_switch(defaultRsw, defaultRcut)
594	>	call setHmatDangerousRcutValue(defaultRcut)
595	>
596	>	haveDefaultCutoffs = .true.
597	>	haveGtypeCutoffMap = .false.
598
599	<	subroutine setCutoffPolicy(cutPolicy)
599	>	end subroutine setCutoffs
600
601	+	subroutine cWasLame()
602	+
603	+	VisitCutoffsAfterComputing = .true.
604	+	return
605	+
606	+	end subroutine cWasLame
607	+
608	+	subroutine setCutoffPolicy(cutPolicy)
609	+
610		integer, intent(in) :: cutPolicy
611	+
612		cutoffPolicy = cutPolicy
613	<	call createGtypeCutoffMap()
613	>	haveCutoffPolicy = .true.
614	>	haveGtypeCutoffMap = .false.
615	>
616		end subroutine setCutoffPolicy
430	–
617
618	<	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()
618	>	subroutine setBoxDipole()
619
620	<	haveSIMvariables = .true.
620	>	do_box_dipole = .true.
621	>
622	>	end subroutine setBoxDipole
623
624	<	return
442	<	end subroutine setSimVariables
624	>	subroutine getBoxDipole( box_dipole )
625
626	+	real(kind=dp), intent(inout), dimension(3) :: box_dipole
627	+
628	+	box_dipole = boxDipole
629	+
630	+	end subroutine getBoxDipole
631	+
632	+	subroutine setElectrostaticMethod( thisESM )
633	+
634	+	integer, intent(in) :: thisESM
635	+
636	+	electrostaticSummationMethod = thisESM
637	+	haveElectrostaticSummationMethod = .true.
638	+
639	+	end subroutine setElectrostaticMethod
640	+
641	+	subroutine setSkinThickness( thisSkin )
642	+
643	+	real(kind=dp), intent(in) :: thisSkin
644	+
645	+	skinThickness = thisSkin
646	+	haveSkinThickness = .true.
647	+	haveGtypeCutoffMap = .false.
648	+
649	+	end subroutine setSkinThickness
650	+
651	+	subroutine setSimVariables()
652	+	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
653	+	SIM_uses_EAM = SimUsesEAM()
654	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
655	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
656	+	SIM_uses_PBC = SimUsesPBC()
657	+	SIM_uses_SC = SimUsesSC()
658	+	SIM_uses_AtomicVirial = SimUsesAtomicVirial()
659	+
660	+	haveSIMvariables = .true.
661	+
662	+	return
663	+	end subroutine setSimVariables
664	+
665		subroutine doReadyCheck(error)
666		integer, intent(out) :: error
446	–
667		integer :: myStatus
668
669		error = 0
670
671		if (.not. haveInteractionHash) then
672	<	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
672	>	call createInteractionHash()
673		endif
674
675		if (.not. haveGtypeCutoffMap) then
676	<	myStatus = 0
463	<	call createGtypeCutoffMap(myStatus)
464	<	if (myStatus .ne. 0) then
465	<	write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
466	<	error = -1
467	<	return
468	<	endif
676	>	call createGtypeCutoffMap()
677		endif
678
679	+	if (VisitCutoffsAfterComputing) then
680	+	call set_switch(largestRcut, largestRcut)
681	+	call setHmatDangerousRcutValue(largestRcut)
682	+	call setCutoffEAM(largestRcut)
683	+	call setCutoffSC(largestRcut)
684	+	VisitCutoffsAfterComputing = .false.
685	+	endif
686	+
687		if (.not. haveSIMvariables) then
688		call setSimVariables()
689		endif
690
475	–	!  if (.not. haveRlist) then
476	–	!     write(default_error, *) 'rList has not been set in doForces!'
477	–	!     error = -1
478	–	!     return
479	–	!  endif
480	–
691		if (.not. haveNeighborList) then
692		write(default_error, *) 'neighbor list has not been initialized in doForces!'
693		error = -1
694		return
695		end if
696	<
696	>
697		if (.not. haveSaneForceField) then
698		write(default_error, *) 'Force Field is not sane in doForces!'
699		error = -1
700		return
701		end if
702	<
702	>
703		#ifdef IS_MPI
704		if (.not. isMPISimSet()) then
705		write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
#	Line 501 | Line 711 | contains
711		end subroutine doReadyCheck
712
713
714	<	subroutine init_FF(thisESM, thisStat)
714	>	subroutine init_FF(thisStat)
715
506	–	integer, intent(in) :: thisESM
507	–	real(kind=dp), intent(in) :: dampingAlpha
716		integer, intent(out) :: thisStat
717		integer :: my_status, nMatches
718		integer, pointer :: MatchList(:) => null()
511	–	real(kind=dp) :: rcut, rrf, rt, dielect
719
720		!! assume things are copacetic, unless they aren't
721		thisStat = 0
722
516	–	electrostaticSummationMethod = thisESM
517	–
723		!! init_FF is called *after* all of the atom types have been
724		!! defined in atype_module using the new_atype subroutine.
725		!!
#	Line 525 | Line 730 | contains
730		FF_uses_Dipoles = .false.
731		FF_uses_GayBerne = .false.
732		FF_uses_EAM = .false.
733	+	FF_uses_SC = .false.
734
735		call getMatchingElementList(atypes, "is_Directional", .true., &
736		nMatches, MatchList)
#	Line 541 | Line 747 | contains
747		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
748		if (nMatches .gt. 0) FF_uses_EAM = .true.
749
750	+	call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
751	+	if (nMatches .gt. 0) FF_uses_SC = .true.
752
753	+
754		haveSaneForceField = .true.
755
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	–
756		if (FF_uses_EAM) then
757		call init_EAM_FF(my_status)
758		if (my_status /= 0) then
#	Line 570 | Line 763 | contains
763		end if
764		endif
765
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	–
766		if (.not. haveNeighborList) then
767		!! Create neighbor lists
768		call expandNeighborList(nLocal, my_status)
#	Line 612 | Line 796 | contains
796
797		!! Stress Tensor
798		real( kind = dp), dimension(9) :: tau
799	<	real ( kind = dp ) :: pot
799	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
800		logical ( kind = 2) :: do_pot_c, do_stress_c
801		logical :: do_pot
802		logical :: do_stress
803		logical :: in_switching_region
804		#ifdef IS_MPI
805	<	real( kind = DP ) :: pot_local
805	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
806		integer :: nAtomsInRow
807		integer :: nAtomsInCol
808		integer :: nprocs
#	Line 631 | Line 815 | contains
815		integer :: istart, iend
816		integer :: ia, jb, atom1, atom2
817		integer :: nlist
818	<	real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
818	>	real( kind = DP ) :: ratmsq, rgrpsq, rgrp, rag, vpair, vij
819		real( kind = DP ) :: sw, dswdr, swderiv, mf
820	<	real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
821	<	real(kind=dp) :: rfpot, mu_i, virial
820	>	real( kind = DP ) :: rVal
821	>	real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij, fg, dag
822	>	real(kind=dp) :: rfpot, mu_i
823	>	real(kind=dp):: rCut
824		integer :: me_i, me_j, n_in_i, n_in_j
825		logical :: is_dp_i
826		integer :: neighborListSize
#	Line 643 | Line 829 | contains
829		integer :: propPack_i, propPack_j
830		integer :: loopStart, loopEnd, loop
831		integer :: iHash
832	<	real(kind=dp) :: listSkin = 1.0
832	>	integer :: i1
833
834	+	!! the variables for the box dipole moment
835	+	#ifdef IS_MPI
836	+	integer :: pChgCount_local
837	+	integer :: nChgCount_local
838	+	real(kind=dp) :: pChg_local
839	+	real(kind=dp) :: nChg_local
840	+	real(kind=dp), dimension(3) :: pChgPos_local
841	+	real(kind=dp), dimension(3) :: nChgPos_local
842	+	real(kind=dp), dimension(3) :: dipVec_local
843	+	#endif
844	+	integer :: pChgCount
845	+	integer :: nChgCount
846	+	real(kind=dp) :: pChg
847	+	real(kind=dp) :: nChg
848	+	real(kind=dp) :: chg_value
849	+	real(kind=dp), dimension(3) :: pChgPos
850	+	real(kind=dp), dimension(3) :: nChgPos
851	+	real(kind=dp), dimension(3) :: dipVec
852	+	real(kind=dp), dimension(3) :: chgVec
853	+
854	+	!! initialize box dipole variables
855	+	if (do_box_dipole) then
856	+	#ifdef IS_MPI
857	+	pChg_local = 0.0_dp
858	+	nChg_local = 0.0_dp
859	+	pChgCount_local = 0
860	+	nChgCount_local = 0
861	+	do i=1, 3
862	+	pChgPos_local = 0.0_dp
863	+	nChgPos_local = 0.0_dp
864	+	dipVec_local = 0.0_dp
865	+	enddo
866	+	#endif
867	+	pChg = 0.0_dp
868	+	nChg = 0.0_dp
869	+	pChgCount = 0
870	+	nChgCount = 0
871	+	chg_value = 0.0_dp
872	+
873	+	do i=1, 3
874	+	pChgPos(i) = 0.0_dp
875	+	nChgPos(i) = 0.0_dp
876	+	dipVec(i) = 0.0_dp
877	+	chgVec(i) = 0.0_dp
878	+	boxDipole(i) = 0.0_dp
879	+	enddo
880	+	endif
881	+
882		!! initialize local variables
883
884		#ifdef IS_MPI
#	Line 707 | Line 941 | contains
941		! (but only on the first time through):
942		if (loop .eq. loopStart) then
943		#ifdef IS_MPI
944	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
944	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
945		update_nlist)
946		#else
947	<	call checkNeighborList(nGroups, q_group, listSkin, &
947	>	call checkNeighborList(nGroups, q_group, skinThickness, &
948		update_nlist)
949		#endif
950		endif
#	Line 768 | Line 1002 | contains
1002		me_j = atid(j)
1003		call get_interatomic_vector(q_group(:,i), &
1004		q_group(:,j), d_grp, rgrpsq)
1005	<	#endif
1005	>	#endif
1006
1007	<	if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
1007	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
1008		if (update_nlist) then
1009		nlist = nlist + 1
1010
#	Line 790 | Line 1024 | contains
1024
1025		list(nlist) = j
1026		endif
1027	+
1028	+	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
1029
1030	<	if (loop .eq. PAIR_LOOP) then
1031	<	vij = 0.0d0
1032	<	fij(1:3) = 0.0d0
1033	<	endif
1034	<
1035	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
1036	<	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:3) = d_grp(1:3)
1052	<	ratmsq = rgrpsq
1053	<	else
1054	<	#ifdef IS_MPI
1055	<	call get_interatomic_vector(q_Row(:,atom1), &
1056	<	q_Col(:,atom2), d_atm, ratmsq)
1030	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
1031	>	if (loop .eq. PAIR_LOOP) then
1032	>	vij = 0.0_dp
1033	>	fij(1) = 0.0_dp
1034	>	fij(2) = 0.0_dp
1035	>	fij(3) = 0.0_dp
1036	>	endif
1037	>
1038	>	call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
1039	>
1040	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
1041	>
1042	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
1043	>
1044	>	atom1 = groupListRow(ia)
1045	>
1046	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1047	>
1048	>	atom2 = groupListCol(jb)
1049	>
1050	>	if (skipThisPair(atom1, atom2))  cycle inner
1051	>
1052	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1053	>	d_atm(1) = d_grp(1)
1054	>	d_atm(2) = d_grp(2)
1055	>	d_atm(3) = d_grp(3)
1056	>	ratmsq = rgrpsq
1057	>	else
1058	>	#ifdef IS_MPI
1059	>	call get_interatomic_vector(q_Row(:,atom1), &
1060	>	q_Col(:,atom2), d_atm, ratmsq)
1061		#else
1062	<	call get_interatomic_vector(q(:,atom1), &
1063	<	q(:,atom2), d_atm, ratmsq)
1062	>	call get_interatomic_vector(q(:,atom1), &
1063	>	q(:,atom2), d_atm, ratmsq)
1064		#endif
1065	<	endif
1066	<
1067	<	if (loop .eq. PREPAIR_LOOP) then
1065	>	endif
1066	>
1067	>	if (loop .eq. PREPAIR_LOOP) then
1068		#ifdef IS_MPI
1069	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1070	<	rgrpsq, d_grp, do_pot, do_stress, &
1071	<	eFrame, A, f, t, pot_local)
1069	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1070	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1071	>	eFrame, A, f, t, pot_local)
1072		#else
1073	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1074	<	rgrpsq, d_grp, do_pot, do_stress, &
1075	<	eFrame, A, f, t, pot)
1073	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1074	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1075	>	eFrame, A, f, t, pot)
1076		#endif
1077	<	else
1077	>	else
1078		#ifdef IS_MPI
1079	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1080	<	do_pot, &
1081	<	eFrame, A, f, t, pot_local, vpair, fpair)
1079	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1080	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1081	>	fpair, d_grp, rgrp, rCut)
1082		#else
1083	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1084	<	do_pot,  &
1085	<	eFrame, A, f, t, pot, vpair, fpair)
1083	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1084	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1085	>	d_grp, rgrp, rCut)
1086		#endif
1087	+	vij = vij + vpair
1088	+	fij(1) = fij(1) + fpair(1)
1089	+	fij(2) = fij(2) + fpair(2)
1090	+	fij(3) = fij(3) + fpair(3)
1091	+	if (do_stress) then
1092	+	call add_stress_tensor(d_atm, fpair, tau)
1093	+	endif
1094	+	endif
1095	+	enddo inner
1096	+	enddo
1097
1098	<	vij = vij + vpair
1099	<	fij(1:3) = fij(1:3) + fpair(1:3)
1100	<	endif
1101	<	enddo inner
1102	<	enddo
1103	<
1104	<	if (loop .eq. PAIR_LOOP) then
1105	<	if (in_switching_region) then
1106	<	swderiv = vij*dswdr/rgrp
1107	<	fij(1) = fij(1) + swderiv*d_grp(1)
1108	<	fij(2) = fij(2) + swderiv*d_grp(2)
1109	<	fij(3) = fij(3) + swderiv*d_grp(3)
1110	<
861	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
862	<	atom1=groupListRow(ia)
863	<	mf = mfactRow(atom1)
1098	>	if (loop .eq. PAIR_LOOP) then
1099	>	if (in_switching_region) then
1100	>	swderiv = vij*dswdr/rgrp
1101	>	fij(1) = fij(1) + swderiv*d_grp(1)
1102	>	fij(2) = fij(2) + swderiv*d_grp(2)
1103	>	fij(3) = fij(3) + swderiv*d_grp(3)
1104	>
1105	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1106	>	atom1=groupListRow(ia)
1107	>	mf = mfactRow(atom1)
1108	>	! fg is the force on atom ia due to cutoff group's
1109	>	! presence in switching region
1110	>	fg = swderiv*d_grp*mf
1111		#ifdef IS_MPI
1112	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1113	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1114	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1112	>	f_Row(1,atom1) = f_Row(1,atom1) + fg(1)
1113	>	f_Row(2,atom1) = f_Row(2,atom1) + fg(2)
1114	>	f_Row(3,atom1) = f_Row(3,atom1) + fg(3)
1115		#else
1116	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1117	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1118	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1116	>	f(1,atom1) = f(1,atom1) + fg(1)
1117	>	f(2,atom1) = f(2,atom1) + fg(2)
1118	>	f(3,atom1) = f(3,atom1) + fg(3)
1119		#endif
1120	<	enddo
1121	<
1122	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1123	<	atom2=groupListCol(jb)
877	<	mf = mfactCol(atom2)
1120	>	if (n_in_i .gt. 1) then
1121	>	if (do_stress.and.SIM_uses_AtomicVirial) then
1122	>	! find the distance between the atom and the center of
1123	>	! the cutoff group:
1124		#ifdef IS_MPI
1125	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1126	<	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
1125	>	call get_interatomic_vector(q_Row(:,atom1), &
1126	>	q_group_Row(:,i), dag, rag)
1127		#else
1128	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1129	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
885	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1128	>	call get_interatomic_vector(q(:,atom1), &
1129	>	q_group(:,i), dag, rag)
1130		#endif
1131	<	enddo
1131	>	call add_stress_tensor(dag,fg,tau)
1132	>	endif
1133	>	endif
1134	>	enddo
1135	>
1136	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1137	>	atom2=groupListCol(jb)
1138	>	mf = mfactCol(atom2)
1139	>	! fg is the force on atom jb due to cutoff group's
1140	>	! presence in switching region
1141	>	fg = -swderiv*d_grp*mf
1142	>	#ifdef IS_MPI
1143	>	f_Col(1,atom2) = f_Col(1,atom2) + fg(1)
1144	>	f_Col(2,atom2) = f_Col(2,atom2) + fg(2)
1145	>	f_Col(3,atom2) = f_Col(3,atom2) + fg(3)
1146	>	#else
1147	>	f(1,atom2) = f(1,atom2) + fg(1)
1148	>	f(2,atom2) = f(2,atom2) + fg(2)
1149	>	f(3,atom2) = f(3,atom2) + fg(3)
1150	>	#endif
1151	>	if (n_in_j .gt. 1) then
1152	>	if (do_stress.and.SIM_uses_AtomicVirial) then
1153	>	! find the distance between the atom and the center of
1154	>	! the cutoff group:
1155	>	#ifdef IS_MPI
1156	>	call get_interatomic_vector(q_Col(:,atom2), &
1157	>	q_group_Col(:,j), dag, rag)
1158	>	#else
1159	>	call get_interatomic_vector(q(:,atom2), &
1160	>	q_group(:,j), dag, rag)
1161	>	#endif
1162	>	call add_stress_tensor(dag,fg,tau)
1163	>	endif
1164	>	endif
1165	>	enddo
1166	>	endif
1167		endif
889	–
890	–	if (do_stress) call add_stress_tensor(d_grp, fij)
1168		endif
1169	<	end if
1169	>	endif
1170		enddo
1171	+
1172		enddo outer
1173
1174		if (update_nlist) then
#	Line 950 | Line 1228 | contains
1228
1229		if (do_pot) then
1230		! scatter/gather pot_row into the members of my column
1231	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1232	<
1231	>	do i = 1,LR_POT_TYPES
1232	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1233	>	end do
1234		! scatter/gather pot_local into all other procs
1235		! add resultant to get total pot
1236		do i = 1, nlocal
1237	<	pot_local = pot_local + pot_Temp(i)
1237	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1238	>	+ pot_Temp(1:LR_POT_TYPES,i)
1239		enddo
1240
1241		pot_Temp = 0.0_DP
1242	<
1243	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1242	>	do i = 1,LR_POT_TYPES
1243	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1244	>	end do
1245		do i = 1, nlocal
1246	<	pot_local = pot_local + pot_Temp(i)
1246	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1247	>	+ pot_Temp(1:LR_POT_TYPES,i)
1248		enddo
1249
1250		endif
1251		#endif
1252
1253	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1253	>	if (SIM_requires_postpair_calc) then
1254	>	do i = 1, nlocal
1255	>
1256	>	! we loop only over the local atoms, so we don't need row and column
1257	>	! lookups for the types
1258	>
1259	>	me_i = atid(i)
1260	>
1261	>	! is the atom electrostatic?  See if it would have an
1262	>	! electrostatic interaction with itself
1263	>	iHash = InteractionHash(me_i,me_i)
1264
1265	<	if (electrostaticSummationMethod == 3) then
974	<
1265	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1266		#ifdef IS_MPI
1267	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1268	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1269	<	do i = 1,nlocal
1270	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1271	<	end do
1267	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1268	>	t, do_pot)
1269	>	#else
1270	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1271	>	t, do_pot)
1272		#endif
1273	<
1274	<	do i = 1, nLocal
1275	<
1276	<	rfpot = 0.0_DP
1273	>	endif
1274	>
1275	>
1276	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1277	>
1278	>	! loop over the excludes to accumulate RF stuff we've
1279	>	! left out of the normal pair loop
1280	>
1281	>	do i1 = 1, nSkipsForAtom(i)
1282	>	j = skipsForAtom(i, i1)
1283	>
1284	>	! prevent overcounting of the skips
1285	>	if (i.lt.j) then
1286	>	call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1287	>	rVal = sqrt(ratmsq)
1288	>	call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1289		#ifdef IS_MPI
1290	<	me_i = atid_row(i)
1290	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1291	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1292		#else
1293	<	me_i = atid(i)
1293	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1294	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1295		#endif
1296	<	iHash = InteractionHash(me_i,me_j)
1297	<
1298	<	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1296	>	endif
1297	>	enddo
1298	>	endif
1299
1300	<	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)
1300	>	if (do_box_dipole) then
1301		#ifdef IS_MPI
1302	<	pot_local = pot_local + rfpot
1302	>	call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1303	>	nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1304	>	pChgCount_local, nChgCount_local)
1305		#else
1306	<	pot = pot + rfpot
1307	<
1306	>	call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1307	>	pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1308		#endif
1309	<	endif
1310	<	enddo
1011	<	endif
1309	>	endif
1310	>	enddo
1311		endif
1312
1014	–
1313		#ifdef IS_MPI
1016	–
1314		if (do_pot) then
1315	<	pot = pot + pot_local
1316	<	!! we assume the c code will do the allreduce to get the total potential
1317	<	!! we could do it right here if we needed to...
1315	>	#ifdef SINGLE_PRECISION
1316	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1317	>	mpi_comm_world,mpi_err)
1318	>	#else
1319	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1320	>	mpi_sum, mpi_comm_world,mpi_err)
1321	>	#endif
1322		endif
1323	+
1324	+	if (do_box_dipole) then
1325
1326	<	if (do_stress) then
1327	<	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1328	<	mpi_comm_world,mpi_err)
1329	<	call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1330	<	mpi_comm_world,mpi_err)
1331	<	endif
1332	<
1326	>	#ifdef SINGLE_PRECISION
1327	>	call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1328	>	mpi_comm_world, mpi_err)
1329	>	call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1330	>	mpi_comm_world, mpi_err)
1331	>	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1332	>	mpi_comm_world, mpi_err)
1333	>	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1334	>	mpi_comm_world, mpi_err)
1335	>	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1336	>	mpi_comm_world, mpi_err)
1337	>	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1338	>	mpi_comm_world, mpi_err)
1339	>	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1340	>	mpi_comm_world, mpi_err)
1341		#else
1342	+	call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1343	+	mpi_comm_world, mpi_err)
1344	+	call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1345	+	mpi_comm_world, mpi_err)
1346	+	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1347	+	mpi_sum, mpi_comm_world, mpi_err)
1348	+	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1349	+	mpi_sum, mpi_comm_world, mpi_err)
1350	+	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1351	+	mpi_sum, mpi_comm_world, mpi_err)
1352	+	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1353	+	mpi_sum, mpi_comm_world, mpi_err)
1354	+	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1355	+	mpi_sum, mpi_comm_world, mpi_err)
1356	+	#endif
1357
1032	–	if (do_stress) then
1033	–	tau = tau_Temp
1034	–	virial = virial_Temp
1358		endif
1359	<
1359	>
1360		#endif
1361
1362	+	if (do_box_dipole) then
1363	+	! first load the accumulated dipole moment (if dipoles were present)
1364	+	boxDipole(1) = dipVec(1)
1365	+	boxDipole(2) = dipVec(2)
1366	+	boxDipole(3) = dipVec(3)
1367	+
1368	+	! now include the dipole moment due to charges
1369	+	! use the lesser of the positive and negative charge totals
1370	+	if (nChg .le. pChg) then
1371	+	chg_value = nChg
1372	+	else
1373	+	chg_value = pChg
1374	+	endif
1375	+
1376	+	! find the average positions
1377	+	if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1378	+	pChgPos = pChgPos / pChgCount
1379	+	nChgPos = nChgPos / nChgCount
1380	+	endif
1381	+
1382	+	! dipole is from the negative to the positive (physics notation)
1383	+	chgVec(1) = pChgPos(1) - nChgPos(1)
1384	+	chgVec(2) = pChgPos(2) - nChgPos(2)
1385	+	chgVec(3) = pChgPos(3) - nChgPos(3)
1386	+
1387	+	boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1388	+	boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1389	+	boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1390	+
1391	+	endif
1392	+
1393		end subroutine do_force_loop
1394
1395		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1396	<	eFrame, A, f, t, pot, vpair, fpair)
1396	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1397
1398	<	real( kind = dp ) :: pot, vpair, sw
1398	>	real( kind = dp ) :: vpair, sw
1399	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1400		real( kind = dp ), dimension(3) :: fpair
1401		real( kind = dp ), dimension(nLocal)   :: mfact
1402		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1052 | Line 1407 | contains
1407		logical, intent(inout) :: do_pot
1408		integer, intent(in) :: i, j
1409		real ( kind = dp ), intent(inout) :: rijsq
1410	<	real ( kind = dp )                :: r
1410	>	real ( kind = dp ), intent(inout) :: r_grp
1411		real ( kind = dp ), intent(inout) :: d(3)
1412	+	real ( kind = dp ), intent(inout) :: d_grp(3)
1413	+	real ( kind = dp ), intent(inout) :: rCut
1414	+	real ( kind = dp ) :: r
1415	+	real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1416		integer :: me_i, me_j
1417	+	integer :: k
1418
1419		integer :: iHash
1420
1421		r = sqrt(rijsq)
1422	<	vpair = 0.0d0
1423	<	fpair(1:3) = 0.0d0
1422	>
1423	>	vpair = 0.0_dp
1424	>	fpair(1:3) = 0.0_dp
1425
1426		#ifdef IS_MPI
1427		me_i = atid_row(i)
#	Line 1071 | Line 1432 | contains
1432		#endif
1433
1434		iHash = InteractionHash(me_i, me_j)
1435	<
1435	>
1436		if ( iand(iHash, LJ_PAIR).ne.0 ) then
1437	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1437	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1438	>	pot(VDW_POT), f, do_pot)
1439		endif
1440	<
1440	>
1441		if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1442	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1443	<	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	<
1442	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1443	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1444		endif
1445	<
1445	>
1446		if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1447		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1448	<	pot, A, f, t, do_pot)
1448	>	pot(HB_POT), A, f, t, do_pot)
1449		endif
1450	<
1450	>
1451		if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1452		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1453	<	pot, A, f, t, do_pot)
1453	>	pot(HB_POT), A, f, t, do_pot)
1454		endif
1455	<
1455	>
1456		if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1457		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1458	<	pot, A, f, t, do_pot)
1458	>	pot(VDW_POT), A, f, t, do_pot)
1459		endif
1460
1461		if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1462	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1463	<	!           pot, A, f, t, do_pot)
1462	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1463	>	pot(VDW_POT), A, f, t, do_pot)
1464		endif
1465	<
1465	>
1466		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1467	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1468	<	do_pot)
1467	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1468	>	pot(METALLIC_POT), f, do_pot)
1469		endif
1470	<
1470	>
1471		if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1472		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1473	<	pot, A, f, t, do_pot)
1473	>	pot(VDW_POT), A, f, t, do_pot)
1474		endif
1475	<
1475	>
1476		if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1477		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1478	<	pot, A, f, t, do_pot)
1478	>	pot(VDW_POT), A, f, t, do_pot)
1479		endif
1480	<
1480	>
1481	>	if ( iand(iHash, SC_PAIR).ne.0 ) then
1482	>	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1483	>	pot(METALLIC_POT), f, do_pot)
1484	>	endif
1485	>
1486		end subroutine do_pair
1487
1488	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1488	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1489		do_pot, do_stress, eFrame, A, f, t, pot)
1490
1491	<	real( kind = dp ) :: pot, sw
1491	>	real( kind = dp ) :: sw
1492	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1493		real( kind = dp ), dimension(9,nLocal) :: eFrame
1494		real (kind=dp), dimension(9,nLocal) :: A
1495		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1137 | Line 1497 | contains
1497
1498		logical, intent(inout) :: do_pot, do_stress
1499		integer, intent(in) :: i, j
1500	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1500	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1501		real ( kind = dp )                :: r, rc
1502		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1503
1504		integer :: me_i, me_j, iHash
1505
1506		r = sqrt(rijsq)
1507	<
1507	>
1508		#ifdef IS_MPI
1509		me_i = atid_row(i)
1510		me_j = atid_col(j)
#	Line 1156 | Line 1516 | contains
1516		iHash = InteractionHash(me_i, me_j)
1517
1518		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1519	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1519	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1520		endif
1521	+
1522	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1523	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1524	+	endif
1525
1526		end subroutine do_prepair
1527
1528
1529		subroutine do_preforce(nlocal,pot)
1530		integer :: nlocal
1531	<	real( kind = dp ) :: pot
1531	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1532
1533		if (FF_uses_EAM .and. SIM_uses_EAM) then
1534	<	call calc_EAM_preforce_Frho(nlocal,pot)
1534	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1535		endif
1536	<
1537	<
1536	>	if (FF_uses_SC .and. SIM_uses_SC) then
1537	>	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1538	>	endif
1539		end subroutine do_preforce
1540
1541
#	Line 1182 | Line 1547 | contains
1547		real( kind = dp ) :: d(3), scaled(3)
1548		integer i
1549
1550	<	d(1:3) = q_j(1:3) - q_i(1:3)
1550	>	d(1) = q_j(1) - q_i(1)
1551	>	d(2) = q_j(2) - q_i(2)
1552	>	d(3) = q_j(3) - q_i(3)
1553
1554		! Wrap back into periodic box if necessary
1555		if ( SIM_uses_PBC ) then
1556
1557		if( .not.boxIsOrthorhombic ) then
1558		! calc the scaled coordinates.
1559	+	! scaled = matmul(HmatInv, d)
1560
1561	<	scaled = matmul(HmatInv, d)
1562	<
1561	>	scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1562	>	scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1563	>	scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1564	>
1565		! wrap the scaled coordinates
1566
1567	<	scaled = scaled  - anint(scaled)
1567	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1568	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1569	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1570
1199	–
1571		! calc the wrapped real coordinates from the wrapped scaled
1572		! coordinates
1573	+	! d = matmul(Hmat,scaled)
1574	+	d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1575	+	d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1576	+	d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1577
1203	–	d = matmul(Hmat,scaled)
1204	–
1578		else
1579		! calc the scaled coordinates.
1580
1581	<	do i = 1, 3
1582	<	scaled(i) = d(i) * HmatInv(i,i)
1583	<
1584	<	! wrap the scaled coordinates
1585	<
1586	<	scaled(i) = scaled(i) - anint(scaled(i))
1581	>	scaled(1) = d(1) * HmatInv(1,1)
1582	>	scaled(2) = d(2) * HmatInv(2,2)
1583	>	scaled(3) = d(3) * HmatInv(3,3)
1584	>
1585	>	! wrap the scaled coordinates
1586	>
1587	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1588	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1589	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1590
1591	<	! calc the wrapped real coordinates from the wrapped scaled
1592	<	! coordinates
1591	>	! calc the wrapped real coordinates from the wrapped scaled
1592	>	! coordinates
1593
1594	<	d(i) = scaled(i)*Hmat(i,i)
1595	<	enddo
1594	>	d(1) = scaled(1)*Hmat(1,1)
1595	>	d(2) = scaled(2)*Hmat(2,2)
1596	>	d(3) = scaled(3)*Hmat(3,3)
1597	>
1598		endif
1599
1600		endif
1601
1602	<	r_sq = dot_product(d,d)
1602	>	r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1603
1604		end subroutine get_interatomic_vector
1605
#	Line 1253 | Line 1631 | contains
1631		pot_Col = 0.0_dp
1632		pot_Temp = 0.0_dp
1633
1256	–	rf_Row = 0.0_dp
1257	–	rf_Col = 0.0_dp
1258	–	rf_Temp = 0.0_dp
1259	–
1634		#endif
1635
1636		if (FF_uses_EAM .and. SIM_uses_EAM) then
1637		call clean_EAM()
1638		endif
1639
1266	–	rf = 0.0_dp
1267	–	tau_Temp = 0.0_dp
1268	–	virial_Temp = 0.0_dp
1640		end subroutine zero_work_arrays
1641
1642		function skipThisPair(atom1, atom2) result(skip_it)
#	Line 1357 | Line 1728 | contains
1728
1729		function FF_RequiresPrepairCalc() result(doesit)
1730		logical :: doesit
1731	<	doesit = FF_uses_EAM
1731	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1732	>	.or. FF_uses_MEAM
1733		end function FF_RequiresPrepairCalc
1734
1363	–	function FF_RequiresPostpairCalc() result(doesit)
1364	–	logical :: doesit
1365	–	if (electrostaticSummationMethod == 3) doesit = .true.
1366	–	end function FF_RequiresPostpairCalc
1367	–
1735		#ifdef PROFILE
1736		function getforcetime() result(totalforcetime)
1737		real(kind=dp) :: totalforcetime
#	Line 1374 | Line 1741 | contains
1741
1742		!! This cleans componets of force arrays belonging only to fortran
1743
1744	<	subroutine add_stress_tensor(dpair, fpair)
1744	>	subroutine add_stress_tensor(dpair, fpair, tau)
1745
1746		real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1747	+	real( kind = dp ), dimension(9), intent(inout) :: tau
1748
1749		! because the d vector is the rj - ri vector, and
1750		! because fx, fy, fz are the force on atom i, we need a
1751		! negative sign here:
1752
1753	<	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1754	<	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1755	<	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1756	<	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1757	<	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1758	<	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1759	<	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1760	<	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1761	<	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1753	>	tau(1) = tau(1) - dpair(1) * fpair(1)
1754	>	tau(2) = tau(2) - dpair(1) * fpair(2)
1755	>	tau(3) = tau(3) - dpair(1) * fpair(3)
1756	>	tau(4) = tau(4) - dpair(2) * fpair(1)
1757	>	tau(5) = tau(5) - dpair(2) * fpair(2)
1758	>	tau(6) = tau(6) - dpair(2) * fpair(3)
1759	>	tau(7) = tau(7) - dpair(3) * fpair(1)
1760	>	tau(8) = tau(8) - dpair(3) * fpair(2)
1761	>	tau(9) = tau(9) - dpair(3) * fpair(3)
1762
1395	–	virial_Temp = virial_Temp + &
1396	–	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1397	–
1763		end subroutine add_stress_tensor
1764
1765		end module doForces

Diff Legend

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