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 2301 by gezelter, Thu Sep 15 22:05:21 2005 UTC vs.
Revision 3132 by chrisfen, Thu May 3 15:52:08 2007 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.44 2005-09-15 22:05:17 gezelter Exp $, $Date: 2005-09-15 22:05:17 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $
48	>	!! @version $Id: doForces.F90,v 1.89 2007-05-03 15:52:08 chrisfen Exp $, $Date: 2007-05-03 15:52:08 $, $Name: not supported by cvs2svn $, $Revision: 1.89 $
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_RF
98	>	logical, save :: FF_uses_SC
99	>	logical, save :: FF_uses_MEAM
100	>
101
102		logical, save :: SIM_uses_DirectionalAtoms
103		logical, save :: SIM_uses_EAM
104	<	logical, save :: SIM_uses_RF
104	>	logical, save :: SIM_uses_SC
105	>	logical, save :: SIM_uses_MEAM
106		logical, save :: SIM_requires_postpair_calc
107		logical, save :: SIM_requires_prepair_calc
108		logical, save :: SIM_uses_PBC
109	+	logical, save :: SIM_uses_AtomicVirial
110
111	<	integer, save :: corrMethod
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
108	–	public :: createInteractionHash
109	–	public :: createGtypeCutoffMap
110	–	public :: getStickyCut
111	–	public :: getStickyPowerCut
112	–	public :: getGayBerneCut
113	–	public :: getEAMCut
114	–	public :: getShapeCut
128
129		#ifdef PROFILE
130		public :: getforcetime
#	Line 124 | 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 134 | 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
139	<	real(kind=dp),save :: rcuti
140	<
155	>	real(kind=dp), dimension(3) :: boxDipole
156	>
157		contains
158
159	<	subroutine createInteractionHash(status)
159	>	subroutine createInteractionHash()
160		integer :: nAtypes
145	–	integer, intent(out) :: status
161		integer :: i
162		integer :: j
163		integer :: iHash
#	Line 154 | 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 161 | 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
166	–	status = 0
167	–
185		if (.not. associated(atypes)) then
186	<	call handleError("atype", "atypes was not present before call of createInteractionHash!")
170	<	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 194 | 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 207 | 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 228 | Line 254 | contains
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)
262		if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
263		if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
#	Line 247 | Line 277 | contains
277		haveInteractionHash = .true.
278		end subroutine createInteractionHash
279
280	<	subroutine createGtypeCutoffMap(stat)
280	>	subroutine createGtypeCutoffMap()
281
252	–	integer, intent(out), optional :: stat
282		logical :: i_is_LJ
283		logical :: i_is_Elect
284		logical :: i_is_Sticky
#	Line 257 | 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
268	–	stat = 0
300		if (.not. haveInteractionHash) then
301	<	call createInteractionHash(myStatus)
271	<	if (myStatus .ne. 0) then
272	<	write(default_error, *) 'createInteractionHash failed in doForces!'
273	<	stat = -1
274	<	return
275	<	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 289 | 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
328	–
329	–	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 359 | 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	<	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
533	<	skin = defaultRlist - defaultRcut
409	<	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
532	>
533	>	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
534
535	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
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)
418	<	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
419	<	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	<	rcuti = 1.0_dp / defaultRcut
426	<	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
434	–
630
631	<	subroutine setSimVariables()
437	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
438	<	SIM_uses_EAM = SimUsesEAM()
439	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
440	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
441	<	SIM_uses_PBC = SimUsesPBC()
442	<	SIM_uses_RF = SimUsesRF()
631	>	subroutine setBoxDipole()
632
633	<	haveSIMvariables = .true.
633	>	do_box_dipole = .true.
634	>
635	>	end subroutine setBoxDipole
636
637	<	return
447	<	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
451	–
680		integer :: myStatus
681
682		error = 0
683
684		if (.not. haveInteractionHash) then
685	<	myStatus = 0
458	<	call createInteractionHash(myStatus)
459	<	if (myStatus .ne. 0) then
460	<	write(default_error, *) 'createInteractionHash failed in doForces!'
461	<	error = -1
462	<	return
463	<	endif
685	>	call createInteractionHash()
686		endif
687
688		if (.not. haveGtypeCutoffMap) then
689	<	myStatus = 0
468	<	call createGtypeCutoffMap(myStatus)
469	<	if (myStatus .ne. 0) then
470	<	write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
471	<	error = -1
472	<	return
473	<	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
480	–	!  if (.not. haveRlist) then
481	–	!     write(default_error, *) 'rList has not been set in doForces!'
482	–	!     error = -1
483	–	!     return
484	–	!  endif
485	–
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 506 | Line 724 | contains
724		end subroutine doReadyCheck
725
726
727	<	subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat)
727	>	subroutine init_FF(thisStat)
728
511	–	logical, intent(in) :: use_RF
512	–	integer, intent(in) :: correctionMethod
513	–	real(kind=dp), intent(in) :: dampingAlpha
729		integer, intent(out) :: thisStat
730		integer :: my_status, nMatches
731		integer, pointer :: MatchList(:) => null()
517	–	real(kind=dp) :: rcut, rrf, rt, dielect
732
733		!! assume things are copacetic, unless they aren't
734		thisStat = 0
735
522	–	!! Fortran's version of a cast:
523	–	FF_uses_RF = use_RF
524	–
525	–
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 533 | 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 549 | 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
555	–	!! check to make sure the FF_uses_RF setting makes sense
556	–
557	–	if (FF_uses_Dipoles) then
558	–	if (FF_uses_RF) then
559	–	dielect = getDielect()
560	–	call initialize_rf(dielect)
561	–	endif
562	–	else
563	–	if ((corrMethod == 3) .or. FF_uses_RF) then
564	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
565	–	thisStat = -1
566	–	haveSaneForceField = .false.
567	–	return
568	–	endif
569	–	endif
570	–
769		if (FF_uses_EAM) then
770		call init_EAM_FF(my_status)
771		if (my_status /= 0) then
#	Line 578 | Line 776 | contains
776		end if
777		endif
778
581	–	if (FF_uses_GayBerne) then
582	–	call check_gb_pair_FF(my_status)
583	–	if (my_status .ne. 0) then
584	–	thisStat = -1
585	–	haveSaneForceField = .false.
586	–	return
587	–	endif
588	–	endif
589	–
779		if (.not. haveNeighborList) then
780		!! Create neighbor lists
781		call expandNeighborList(nLocal, my_status)
#	Line 620 | 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 639 | 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 651 | 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 715 | 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 776 | 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 798 | 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)
885	<	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
889	<	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
893	<	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
897	–
898	–	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 958 | Line 1247 | contains
1247
1248		if (do_pot) then
1249		! scatter/gather pot_row into the members of my column
1250	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1251	<
1250	>	do i = 1,LR_POT_TYPES
1251	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1252	>	end do
1253		! scatter/gather pot_local into all other procs
1254		! add resultant to get total pot
1255		do i = 1, nlocal
1256	<	pot_local = pot_local + pot_Temp(i)
1256	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1257	>	+ pot_Temp(1:LR_POT_TYPES,i)
1258		enddo
1259
1260		pot_Temp = 0.0_DP
1261	<
1262	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1261	>	do i = 1,LR_POT_TYPES
1262	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1263	>	end do
1264		do i = 1, nlocal
1265	<	pot_local = pot_local + pot_Temp(i)
1265	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1266	>	+ pot_Temp(1:LR_POT_TYPES,i)
1267		enddo
1268
1269		endif
1270		#endif
1271
1272	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1272	>	if (SIM_requires_postpair_calc) then
1273	>	do i = 1, nlocal
1274	>
1275	>	! we loop only over the local atoms, so we don't need row and column
1276	>	! lookups for the types
1277	>
1278	>	me_i = atid(i)
1279	>
1280	>	! is the atom electrostatic?  See if it would have an
1281	>	! electrostatic interaction with itself
1282	>	iHash = InteractionHash(me_i,me_i)
1283
1284	<	if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
982	<
1284	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1285		#ifdef IS_MPI
1286	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1287	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1288	<	do i = 1,nlocal
1289	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1290	<	end do
1286	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1287	>	t, do_pot)
1288	>	#else
1289	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1290	>	t, do_pot)
1291		#endif
1292	<
1293	<	do i = 1, nLocal
1294	<
1295	<	rfpot = 0.0_DP
1292	>	endif
1293	>
1294	>
1295	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1296	>
1297	>	! loop over the excludes to accumulate RF stuff we've
1298	>	! left out of the normal pair loop
1299	>
1300	>	do i1 = 1, nSkipsForAtom(i)
1301	>	j = skipsForAtom(i, i1)
1302	>
1303	>	! prevent overcounting of the skips
1304	>	if (i.lt.j) then
1305	>	call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1306	>	rVal = sqrt(ratmsq)
1307	>	call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1308		#ifdef IS_MPI
1309	<	me_i = atid_row(i)
1309	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1310	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1311		#else
1312	<	me_i = atid(i)
1312	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1313	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1314		#endif
1315	<	iHash = InteractionHash(me_i,me_j)
1316	<
1317	<	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1315	>	endif
1316	>	enddo
1317	>	endif
1318
1319	<	mu_i = getDipoleMoment(me_i)
1004	<
1005	<	!! The reaction field needs to include a self contribution
1006	<	!! to the field:
1007	<	call accumulate_self_rf(i, mu_i, eFrame)
1008	<	!! Get the reaction field contribution to the
1009	<	!! potential and torques:
1010	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1319	>	if (do_box_dipole) then
1320		#ifdef IS_MPI
1321	<	pot_local = pot_local + rfpot
1321	>	call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1322	>	nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1323	>	pChgCount_local, nChgCount_local)
1324		#else
1325	<	pot = pot + rfpot
1326	<
1325	>	call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1326	>	pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1327		#endif
1328	<	endif
1329	<	enddo
1019	<	endif
1328	>	endif
1329	>	enddo
1330		endif
1331
1022	–
1332		#ifdef IS_MPI
1024	–
1333		if (do_pot) then
1334	<	pot = pot + pot_local
1335	<	!! we assume the c code will do the allreduce to get the total potential
1336	<	!! we could do it right here if we needed to...
1334	>	#ifdef SINGLE_PRECISION
1335	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1336	>	mpi_comm_world,mpi_err)
1337	>	#else
1338	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1339	>	mpi_sum, mpi_comm_world,mpi_err)
1340	>	#endif
1341		endif
1342	+
1343	+	if (do_box_dipole) then
1344
1345	<	if (do_stress) then
1346	<	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1347	<	mpi_comm_world,mpi_err)
1348	<	call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1349	<	mpi_comm_world,mpi_err)
1350	<	endif
1351	<
1345	>	#ifdef SINGLE_PRECISION
1346	>	call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1347	>	mpi_comm_world, mpi_err)
1348	>	call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1349	>	mpi_comm_world, mpi_err)
1350	>	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1351	>	mpi_comm_world, mpi_err)
1352	>	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1353	>	mpi_comm_world, mpi_err)
1354	>	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1355	>	mpi_comm_world, mpi_err)
1356	>	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1357	>	mpi_comm_world, mpi_err)
1358	>	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1359	>	mpi_comm_world, mpi_err)
1360		#else
1361	+	call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1362	+	mpi_comm_world, mpi_err)
1363	+	call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1364	+	mpi_comm_world, mpi_err)
1365	+	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1366	+	mpi_sum, mpi_comm_world, mpi_err)
1367	+	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1368	+	mpi_sum, mpi_comm_world, mpi_err)
1369	+	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1370	+	mpi_sum, mpi_comm_world, mpi_err)
1371	+	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1372	+	mpi_sum, mpi_comm_world, mpi_err)
1373	+	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1374	+	mpi_sum, mpi_comm_world, mpi_err)
1375	+	#endif
1376
1040	–	if (do_stress) then
1041	–	tau = tau_Temp
1042	–	virial = virial_Temp
1377		endif
1378	<
1378	>
1379		#endif
1380
1381	+	if (do_box_dipole) then
1382	+	! first load the accumulated dipole moment (if dipoles were present)
1383	+	boxDipole(1) = dipVec(1)
1384	+	boxDipole(2) = dipVec(2)
1385	+	boxDipole(3) = dipVec(3)
1386	+
1387	+	! now include the dipole moment due to charges
1388	+	! use the lesser of the positive and negative charge totals
1389	+	if (nChg .le. pChg) then
1390	+	chg_value = nChg
1391	+	else
1392	+	chg_value = pChg
1393	+	endif
1394	+
1395	+	! find the average positions
1396	+	if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1397	+	pChgPos = pChgPos / pChgCount
1398	+	nChgPos = nChgPos / nChgCount
1399	+	endif
1400	+
1401	+	! dipole is from the negative to the positive (physics notation)
1402	+	chgVec(1) = pChgPos(1) - nChgPos(1)
1403	+	chgVec(2) = pChgPos(2) - nChgPos(2)
1404	+	chgVec(3) = pChgPos(3) - nChgPos(3)
1405	+
1406	+	boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1407	+	boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1408	+	boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1409	+
1410	+	endif
1411	+
1412		end subroutine do_force_loop
1413
1414		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1415	<	eFrame, A, f, t, pot, vpair, fpair)
1415	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1416
1417	<	real( kind = dp ) :: pot, vpair, sw
1417	>	real( kind = dp ) :: vpair, sw
1418	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1419		real( kind = dp ), dimension(3) :: fpair
1420		real( kind = dp ), dimension(nLocal)   :: mfact
1421		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1060 | Line 1426 | contains
1426		logical, intent(inout) :: do_pot
1427		integer, intent(in) :: i, j
1428		real ( kind = dp ), intent(inout) :: rijsq
1429	<	real ( kind = dp )                :: r
1429	>	real ( kind = dp ), intent(inout) :: r_grp
1430		real ( kind = dp ), intent(inout) :: d(3)
1431	+	real ( kind = dp ), intent(inout) :: d_grp(3)
1432	+	real ( kind = dp ), intent(inout) :: rCut
1433	+	real ( kind = dp ) :: r
1434	+	real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1435		integer :: me_i, me_j
1436	+	integer :: k
1437
1438		integer :: iHash
1439
1440		r = sqrt(rijsq)
1441	<	vpair = 0.0d0
1442	<	fpair(1:3) = 0.0d0
1441	>
1442	>	vpair = 0.0_dp
1443	>	fpair(1:3) = 0.0_dp
1444
1445		#ifdef IS_MPI
1446		me_i = atid_row(i)
#	Line 1079 | Line 1451 | contains
1451		#endif
1452
1453		iHash = InteractionHash(me_i, me_j)
1454	<
1454	>
1455		if ( iand(iHash, LJ_PAIR).ne.0 ) then
1456	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1456	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1457	>	pot(VDW_POT), f, do_pot)
1458		endif
1459	<
1459	>
1460		if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1461	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1462	<	pot, eFrame, f, t, do_pot, corrMethod, rcuti)
1090	<
1091	<	if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
1092	<
1093	<	! CHECK ME (RF needs to know about all electrostatic types)
1094	<	call accumulate_rf(i, j, r, eFrame, sw)
1095	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1096	<	endif
1097	<
1461	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1462	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1463		endif
1464	<
1464	>
1465		if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1466		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1467	<	pot, A, f, t, do_pot)
1467	>	pot(HB_POT), A, f, t, do_pot)
1468		endif
1469	<
1469	>
1470		if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1471		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1472	<	pot, A, f, t, do_pot)
1472	>	pot(HB_POT), A, f, t, do_pot)
1473		endif
1474	<
1474	>
1475		if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1476		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1477	<	pot, A, f, t, do_pot)
1477	>	pot(VDW_POT), A, f, t, do_pot)
1478		endif
1479
1480		if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1481	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1482	<	!           pot, A, f, t, do_pot)
1481	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1482	>	pot(VDW_POT), A, f, t, do_pot)
1483		endif
1484	<
1484	>
1485		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1486	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1487	<	do_pot)
1486	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1487	>	pot(METALLIC_POT), f, do_pot)
1488		endif
1489	<
1489	>
1490		if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1491		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1492	<	pot, A, f, t, do_pot)
1492	>	pot(VDW_POT), A, f, t, do_pot)
1493		endif
1494	<
1494	>
1495		if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1496		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1497	<	pot, A, f, t, do_pot)
1497	>	pot(VDW_POT), A, f, t, do_pot)
1498		endif
1499	<
1499	>
1500	>	if ( iand(iHash, SC_PAIR).ne.0 ) then
1501	>	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1502	>	pot(METALLIC_POT), f, do_pot)
1503	>	endif
1504	>
1505		end subroutine do_pair
1506
1507	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1507	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1508		do_pot, do_stress, eFrame, A, f, t, pot)
1509
1510	<	real( kind = dp ) :: pot, sw
1510	>	real( kind = dp ) :: sw
1511	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1512		real( kind = dp ), dimension(9,nLocal) :: eFrame
1513		real (kind=dp), dimension(9,nLocal) :: A
1514		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1145 | Line 1516 | contains
1516
1517		logical, intent(inout) :: do_pot, do_stress
1518		integer, intent(in) :: i, j
1519	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1519	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1520		real ( kind = dp )                :: r, rc
1521		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1522
1523		integer :: me_i, me_j, iHash
1524
1525		r = sqrt(rijsq)
1526	<
1526	>
1527		#ifdef IS_MPI
1528		me_i = atid_row(i)
1529		me_j = atid_col(j)
#	Line 1164 | Line 1535 | contains
1535		iHash = InteractionHash(me_i, me_j)
1536
1537		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1538	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1538	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1539		endif
1540	+
1541	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1542	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1543	+	endif
1544
1545		end subroutine do_prepair
1546
1547
1548		subroutine do_preforce(nlocal,pot)
1549		integer :: nlocal
1550	<	real( kind = dp ) :: pot
1550	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1551
1552		if (FF_uses_EAM .and. SIM_uses_EAM) then
1553	<	call calc_EAM_preforce_Frho(nlocal,pot)
1553	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1554		endif
1555	<
1556	<
1555	>	if (FF_uses_SC .and. SIM_uses_SC) then
1556	>	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1557	>	endif
1558		end subroutine do_preforce
1559
1560
#	Line 1190 | Line 1566 | contains
1566		real( kind = dp ) :: d(3), scaled(3)
1567		integer i
1568
1569	<	d(1:3) = q_j(1:3) - q_i(1:3)
1569	>	d(1) = q_j(1) - q_i(1)
1570	>	d(2) = q_j(2) - q_i(2)
1571	>	d(3) = q_j(3) - q_i(3)
1572
1573		! Wrap back into periodic box if necessary
1574		if ( SIM_uses_PBC ) then
1575
1576		if( .not.boxIsOrthorhombic ) then
1577		! calc the scaled coordinates.
1578	+	! scaled = matmul(HmatInv, d)
1579
1580	<	scaled = matmul(HmatInv, d)
1581	<
1580	>	scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1581	>	scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1582	>	scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1583	>
1584		! wrap the scaled coordinates
1585
1586	<	scaled = scaled  - anint(scaled)
1586	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1587	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1588	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1589
1207	–
1590		! calc the wrapped real coordinates from the wrapped scaled
1591		! coordinates
1592	+	! d = matmul(Hmat,scaled)
1593	+	d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1594	+	d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1595	+	d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1596
1211	–	d = matmul(Hmat,scaled)
1212	–
1597		else
1598		! calc the scaled coordinates.
1599
1600	<	do i = 1, 3
1601	<	scaled(i) = d(i) * HmatInv(i,i)
1602	<
1603	<	! wrap the scaled coordinates
1600	>	scaled(1) = d(1) * HmatInv(1,1)
1601	>	scaled(2) = d(2) * HmatInv(2,2)
1602	>	scaled(3) = d(3) * HmatInv(3,3)
1603	>
1604	>	! wrap the scaled coordinates
1605	>
1606	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1607	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1608	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1609
1610	<	scaled(i) = scaled(i) - anint(scaled(i))
1610	>	! calc the wrapped real coordinates from the wrapped scaled
1611	>	! coordinates
1612
1613	<	! calc the wrapped real coordinates from the wrapped scaled
1614	<	! coordinates
1613	>	d(1) = scaled(1)*Hmat(1,1)
1614	>	d(2) = scaled(2)*Hmat(2,2)
1615	>	d(3) = scaled(3)*Hmat(3,3)
1616
1226	–	d(i) = scaled(i)*Hmat(i,i)
1227	–	enddo
1617		endif
1618
1619		endif
1620
1621	<	r_sq = dot_product(d,d)
1621	>	r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1622
1623		end subroutine get_interatomic_vector
1624
#	Line 1261 | Line 1650 | contains
1650		pot_Col = 0.0_dp
1651		pot_Temp = 0.0_dp
1652
1264	–	rf_Row = 0.0_dp
1265	–	rf_Col = 0.0_dp
1266	–	rf_Temp = 0.0_dp
1267	–
1653		#endif
1654
1655		if (FF_uses_EAM .and. SIM_uses_EAM) then
1656		call clean_EAM()
1657		endif
1658
1274	–	rf = 0.0_dp
1275	–	tau_Temp = 0.0_dp
1276	–	virial_Temp = 0.0_dp
1659		end subroutine zero_work_arrays
1660
1661		function skipThisPair(atom1, atom2) result(skip_it)
#	Line 1365 | Line 1747 | contains
1747
1748		function FF_RequiresPrepairCalc() result(doesit)
1749		logical :: doesit
1750	<	doesit = FF_uses_EAM
1750	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1751	>	.or. FF_uses_MEAM
1752		end function FF_RequiresPrepairCalc
1753
1371	–	function FF_RequiresPostpairCalc() result(doesit)
1372	–	logical :: doesit
1373	–	doesit = FF_uses_RF
1374	–	if (corrMethod == 3) doesit = .true.
1375	–	end function FF_RequiresPostpairCalc
1376	–
1754		#ifdef PROFILE
1755		function getforcetime() result(totalforcetime)
1756		real(kind=dp) :: totalforcetime
#	Line 1383 | Line 1760 | contains
1760
1761		!! This cleans componets of force arrays belonging only to fortran
1762
1763	<	subroutine add_stress_tensor(dpair, fpair)
1763	>	subroutine add_stress_tensor(dpair, fpair, tau)
1764
1765		real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1766	+	real( kind = dp ), dimension(9), intent(inout) :: tau
1767
1768		! because the d vector is the rj - ri vector, and
1769		! because fx, fy, fz are the force on atom i, we need a
1770		! negative sign here:
1771
1772	<	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1773	<	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1774	<	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1775	<	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1776	<	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1777	<	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1778	<	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1779	<	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1780	<	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1772	>	tau(1) = tau(1) - dpair(1) * fpair(1)
1773	>	tau(2) = tau(2) - dpair(1) * fpair(2)
1774	>	tau(3) = tau(3) - dpair(1) * fpair(3)
1775	>	tau(4) = tau(4) - dpair(2) * fpair(1)
1776	>	tau(5) = tau(5) - dpair(2) * fpair(2)
1777	>	tau(6) = tau(6) - dpair(2) * fpair(3)
1778	>	tau(7) = tau(7) - dpair(3) * fpair(1)
1779	>	tau(8) = tau(8) - dpair(3) * fpair(2)
1780	>	tau(9) = tau(9) - dpair(3) * fpair(3)
1781
1404	–	virial_Temp = virial_Temp + &
1405	–	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1406	–
1782		end subroutine add_stress_tensor
1783
1784		end module doForces

Diff Legend

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