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 2266 by chuckv, Thu Jul 28 22:12:45 2005 UTC vs.
Revision 2592 by gezelter, Thu Feb 16 21:40:20 2006 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.24 2005-07-28 22:12:45 chuckv Exp $, $Date: 2005-07-28 22:12:45 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
48	>	!! @version $Id: doForces.F90,v 1.76 2006-02-16 21:40:20 gezelter Exp $, $Date: 2006-02-16 21:40:20 $, $Name: not supported by cvs2svn $, $Revision: 1.76 $
49
50
51		module doForces
#	Line 58 | Line 58 | module doForces
58		use lj
59		use sticky
60		use electrostatic_module
61	<	use reaction_field
62	<	use gb_pair
61	>	use gayberne
62		use shapes
63		use vector_class
64		use eam
65	+	use suttonchen
66		use status
67		#ifdef IS_MPI
68		use mpiSimulation
#	Line 73 | Line 73 | module doForces
73
74		#define __FORTRAN90
75		#include "UseTheForce/fSwitchingFunction.h"
76	+	#include "UseTheForce/fCutoffPolicy.h"
77		#include "UseTheForce/DarkSide/fInteractionMap.h"
78	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79
80	+
81		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
82		INTEGER, PARAMETER:: PAIR_LOOP    = 2
83
81	–	logical, save :: haveRlist = .false.
84		logical, save :: haveNeighborList = .false.
85		logical, save :: haveSIMvariables = .false.
86		logical, save :: haveSaneForceField = .false.
87	<	logical, save :: haveInteractionMap = .false.
87	>	logical, save :: haveInteractionHash = .false.
88	>	logical, save :: haveGtypeCutoffMap = .false.
89	>	logical, save :: haveDefaultCutoffs = .false.
90	>	logical, save :: haveSkinThickness = .false.
91	>	logical, save :: haveElectrostaticSummationMethod = .false.
92	>	logical, save :: haveCutoffPolicy = .false.
93	>	logical, save :: VisitCutoffsAfterComputing = .false.
94
95		logical, save :: FF_uses_DirectionalAtoms
88	–	logical, save :: FF_uses_LennardJones
89	–	logical, save :: FF_uses_Electrostatics
90	–	logical, save :: FF_uses_Charges
96		logical, save :: FF_uses_Dipoles
92	–	logical, save :: FF_uses_Quadrupoles
93	–	logical, save :: FF_uses_Sticky
94	–	logical, save :: FF_uses_StickyPower
97		logical, save :: FF_uses_GayBerne
98		logical, save :: FF_uses_EAM
99	<	logical, save :: FF_uses_Shapes
100	<	logical, save :: FF_uses_FLARB
101	<	logical, save :: FF_uses_RF
99	>	logical, save :: FF_uses_SC
100	>	logical, save :: FF_uses_MEAM
101	>
102
103		logical, save :: SIM_uses_DirectionalAtoms
102	–	logical, save :: SIM_uses_LennardJones
103	–	logical, save :: SIM_uses_Electrostatics
104	–	logical, save :: SIM_uses_Charges
105	–	logical, save :: SIM_uses_Dipoles
106	–	logical, save :: SIM_uses_Quadrupoles
107	–	logical, save :: SIM_uses_Sticky
108	–	logical, save :: SIM_uses_StickyPower
109	–	logical, save :: SIM_uses_GayBerne
104		logical, save :: SIM_uses_EAM
105	<	logical, save :: SIM_uses_Shapes
106	<	logical, save :: SIM_uses_FLARB
113	<	logical, save :: SIM_uses_RF
105	>	logical, save :: SIM_uses_SC
106	>	logical, save :: SIM_uses_MEAM
107		logical, save :: SIM_requires_postpair_calc
108		logical, save :: SIM_requires_prepair_calc
109		logical, save :: SIM_uses_PBC
117	–	logical, save :: SIM_uses_molecular_cutoffs
110
111	<	!!!GO AWAY---------
112	<	!!!!!real(kind=dp), save :: rlist, rlistsq
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 :: setCutoffs
120	+	public :: cWasLame
121	+	public :: setElectrostaticMethod
122	+	public :: setCutoffPolicy
123	+	public :: setSkinThickness
124		public :: do_force_loop
124	–	!  public :: setRlistDF
125	–	!public :: addInteraction
126	–	!public :: setInteractionHash
127	–	!public :: getInteractionHash
128	–	public :: createInteractionMap
129	–	public :: createRcuts
125
126		#ifdef PROFILE
127		public :: getforcetime
#	Line 134 | Line 129 | module doForces
129		real :: forceTimeInitial, forceTimeFinal
130		integer :: nLoops
131		#endif
137	–
138	–	type, public :: Interaction
139	–	integer :: InteractionHash
140	–	real(kind=dp) :: rList = 0.0_dp
141	–	real(kind=dp) :: rListSq = 0.0_dp
142	–	end type Interaction
132
133	<	type(Interaction), dimension(:,:),allocatable :: InteractionMap
134	<
133	>	!! Variables for cutoff mapping and interaction mapping
134	>	! Bit hash to determine pair-pair interactions.
135	>	integer, dimension(:,:), allocatable :: InteractionHash
136	>	real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
137	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
138	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
139
140	<
140	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
141	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
142	>
143	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
144	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
145	>	type ::gtypeCutoffs
146	>	real(kind=dp) :: rcut
147	>	real(kind=dp) :: rcutsq
148	>	real(kind=dp) :: rlistsq
149	>	end type gtypeCutoffs
150	>	type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
151	>
152		contains
153
154	<
151	<	subroutine createInteractionMap(status)
154	>	subroutine createInteractionHash()
155		integer :: nAtypes
153	–	integer, intent(out) :: status
156		integer :: i
157		integer :: j
158	<	integer :: ihash
159	<	real(kind=dp) :: myRcut
158	<	! Test Types
158	>	integer :: iHash
159	>	!! Test Types
160		logical :: i_is_LJ
161		logical :: i_is_Elect
162		logical :: i_is_Sticky
#	Line 163 | Line 164 | contains
164		logical :: i_is_GB
165		logical :: i_is_EAM
166		logical :: i_is_Shape
167	+	logical :: i_is_SC
168	+	logical :: i_is_MEAM
169		logical :: j_is_LJ
170		logical :: j_is_Elect
171		logical :: j_is_Sticky
#	Line 170 | Line 173 | contains
173		logical :: j_is_GB
174		logical :: j_is_EAM
175		logical :: j_is_Shape
176	<
177	<	status = 0
178	<
176	>	logical :: j_is_SC
177	>	logical :: j_is_MEAM
178	>	real(kind=dp) :: myRcut
179	>
180		if (.not. associated(atypes)) then
181	<	call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
178	<	status = -1
181	>	call handleError("doForces", "atypes was not present before call of createInteractionHash!")
182		return
183		endif
184
185		nAtypes = getSize(atypes)
186
187		if (nAtypes == 0) then
188	<	status = -1
188	>	call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
189		return
190		end if
191
192	<	if (.not. allocated(InteractionMap)) then
193	<	allocate(InteractionMap(nAtypes,nAtypes))
192	>	if (.not. allocated(InteractionHash)) then
193	>	allocate(InteractionHash(nAtypes,nAtypes))
194	>	else
195	>	deallocate(InteractionHash)
196	>	allocate(InteractionHash(nAtypes,nAtypes))
197		endif
198	+
199	+	if (.not. allocated(atypeMaxCutoff)) then
200	+	allocate(atypeMaxCutoff(nAtypes))
201	+	else
202	+	deallocate(atypeMaxCutoff)
203	+	allocate(atypeMaxCutoff(nAtypes))
204	+	endif
205
206		do i = 1, nAtypes
207		call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
#	Line 198 | Line 211 | contains
211		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
212		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
213		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
214	+	call getElementProperty(atypes, i, "is_SC", i_is_SC)
215	+	call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
216
217		do j = i, nAtypes
218
#	Line 211 | Line 226 | contains
226		call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
227		call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
228		call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
229	+	call getElementProperty(atypes, j, "is_SC", j_is_SC)
230	+	call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
231
232		if (i_is_LJ .and. j_is_LJ) then
233		iHash = ior(iHash, LJ_PAIR)
#	Line 232 | Line 249 | contains
249		iHash = ior(iHash, EAM_PAIR)
250		endif
251
252	+	if (i_is_SC .and. j_is_SC) then
253	+	iHash = ior(iHash, SC_PAIR)
254	+	endif
255	+
256		if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
257		if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
258		if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
#	Line 241 | Line 262 | contains
262		if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
263
264
265	<	InteractionMap(i,j)%InteractionHash = iHash
266	<	InteractionMap(j,i)%InteractionHash = iHash
265	>	InteractionHash(i,j) = iHash
266	>	InteractionHash(j,i) = iHash
267
268		end do
269
270		end do
250	–	end subroutine createInteractionMap
271
272	<	! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff.
273	<	subroutine createRcuts(defaultRList,stat)
254	<	real(kind=dp), intent(in), optional :: defaultRList
255	<	integer :: iMap
256	<	integer :: map_i,map_j
257	<	real(kind=dp) :: thisRCut = 0.0_dp
258	<	real(kind=dp) :: actualCutoff = 0.0_dp
259	<	integer, intent(out) :: stat
260	<	integer :: nAtypes
261	<	integer :: myStatus
272	>	haveInteractionHash = .true.
273	>	end subroutine createInteractionHash
274
275	<	stat = 0
264	<	if (.not. haveInteractionMap) then
275	>	subroutine createGtypeCutoffMap()
276
277	<	call createInteractionMap(myStatus)
278	<
279	<	if (myStatus .ne. 0) then
280	<	write(default_error, *) 'createInteractionMap failed in doForces!'
281	<	stat = -1
282	<	return
277	>	logical :: i_is_LJ
278	>	logical :: i_is_Elect
279	>	logical :: i_is_Sticky
280	>	logical :: i_is_StickyP
281	>	logical :: i_is_GB
282	>	logical :: i_is_EAM
283	>	logical :: i_is_Shape
284	>	logical :: i_is_SC
285	>	logical :: GtypeFound
286	>
287	>	integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
288	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
289	>	integer :: nGroupsInRow
290	>	integer :: nGroupsInCol
291	>	integer :: nGroupTypesRow,nGroupTypesCol
292	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
293	>	real(kind=dp) :: biggestAtypeCutoff
294	>
295	>	if (.not. haveInteractionHash) then
296	>	call createInteractionHash()
297	>	endif
298	>	#ifdef IS_MPI
299	>	nGroupsInRow = getNgroupsInRow(plan_group_row)
300	>	nGroupsInCol = getNgroupsInCol(plan_group_col)
301	>	#endif
302	>	nAtypes = getSize(atypes)
303	>	! Set all of the initial cutoffs to zero.
304	>	atypeMaxCutoff = 0.0_dp
305	>	do i = 1, nAtypes
306	>	if (SimHasAtype(i)) then
307	>	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
308	>	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
309	>	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
310	>	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
311	>	call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
312	>	call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
313	>	call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
314	>	call getElementProperty(atypes, i, "is_SC", i_is_SC)
315	>
316	>	if (haveDefaultCutoffs) then
317	>	atypeMaxCutoff(i) = defaultRcut
318	>	else
319	>	if (i_is_LJ) then
320	>	thisRcut = getSigma(i) * 2.5_dp
321	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322	>	endif
323	>	if (i_is_Elect) then
324	>	thisRcut = defaultRcut
325	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326	>	endif
327	>	if (i_is_Sticky) then
328	>	thisRcut = getStickyCut(i)
329	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330	>	endif
331	>	if (i_is_StickyP) then
332	>	thisRcut = getStickyPowerCut(i)
333	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334	>	endif
335	>	if (i_is_GB) then
336	>	thisRcut = getGayBerneCut(i)
337	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
338	>	endif
339	>	if (i_is_EAM) then
340	>	thisRcut = getEAMCut(i)
341	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
342	>	endif
343	>	if (i_is_Shape) then
344	>	thisRcut = getShapeCut(i)
345	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
346	>	endif
347	>	if (i_is_SC) then
348	>	thisRcut = getSCCut(i)
349	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
350	>	endif
351	>	endif
352	>
353	>	if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
354	>	biggestAtypeCutoff = atypeMaxCutoff(i)
355	>	endif
356	>
357		endif
358	+	enddo
359	+
360	+	istart = 1
361	+	jstart = 1
362	+	#ifdef IS_MPI
363	+	iend = nGroupsInRow
364	+	jend = nGroupsInCol
365	+	#else
366	+	iend = nGroups
367	+	jend = nGroups
368	+	#endif
369	+
370	+	!! allocate the groupToGtype and gtypeMaxCutoff here.
371	+	if(.not.allocated(groupToGtypeRow)) then
372	+	!  allocate(groupToGtype(iend))
373	+	allocate(groupToGtypeRow(iend))
374	+	else
375	+	deallocate(groupToGtypeRow)
376	+	allocate(groupToGtypeRow(iend))
377		endif
378	+	if(.not.allocated(groupMaxCutoffRow)) then
379	+	allocate(groupMaxCutoffRow(iend))
380	+	else
381	+	deallocate(groupMaxCutoffRow)
382	+	allocate(groupMaxCutoffRow(iend))
383	+	end if
384
385	+	if(.not.allocated(gtypeMaxCutoffRow)) then
386	+	allocate(gtypeMaxCutoffRow(iend))
387	+	else
388	+	deallocate(gtypeMaxCutoffRow)
389	+	allocate(gtypeMaxCutoffRow(iend))
390	+	endif
391
392	<	nAtypes = getSize(atypes)
393	<	! If we pass a default rcut, set all atypes to that cutoff distance
394	<	if(present(defaultRList)) then
395	<	InteractionMap(:,:)%rList = defaultRList
396	<	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
397	<	haveRlist = .true.
398	<	return
392	>
393	>	#ifdef IS_MPI
394	>	! We only allocate new storage if we are in MPI because Ncol /= Nrow
395	>	if(.not.associated(groupToGtypeCol)) then
396	>	allocate(groupToGtypeCol(jend))
397	>	else
398	>	deallocate(groupToGtypeCol)
399	>	allocate(groupToGtypeCol(jend))
400		end if
401
402	<	do map_i = 1,nAtypes
403	<	do map_j = map_i,nAtypes
404	<	iMap = InteractionMap(map_i, map_j)%InteractionHash
405	<
406	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
407	<	!            thisRCut = getLJCutOff(map_i,map_j)
408	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
409	<	endif
410	<
411	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
412	<	!            thisRCut = getElectrostaticCutOff(map_i,map_j)
413	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
402	>	if(.not.associated(groupMaxCutoffCol)) then
403	>	allocate(groupMaxCutoffCol(jend))
404	>	else
405	>	deallocate(groupMaxCutoffCol)
406	>	allocate(groupMaxCutoffCol(jend))
407	>	end if
408	>	if(.not.associated(gtypeMaxCutoffCol)) then
409	>	allocate(gtypeMaxCutoffCol(jend))
410	>	else
411	>	deallocate(gtypeMaxCutoffCol)
412	>	allocate(gtypeMaxCutoffCol(jend))
413	>	end if
414	>
415	>	groupMaxCutoffCol = 0.0_dp
416	>	gtypeMaxCutoffCol = 0.0_dp
417	>
418	>	#endif
419	>	groupMaxCutoffRow = 0.0_dp
420	>	gtypeMaxCutoffRow = 0.0_dp
421	>
422	>
423	>	!! first we do a single loop over the cutoff groups to find the
424	>	!! largest cutoff for any atypes present in this group.  We also
425	>	!! create gtypes at this point.
426	>
427	>	tol = 1.0d-6
428	>	nGroupTypesRow = 0
429	>	nGroupTypesCol = 0
430	>	do i = istart, iend
431	>	n_in_i = groupStartRow(i+1) - groupStartRow(i)
432	>	groupMaxCutoffRow(i) = 0.0_dp
433	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
434	>	atom1 = groupListRow(ia)
435	>	#ifdef IS_MPI
436	>	me_i = atid_row(atom1)
437	>	#else
438	>	me_i = atid(atom1)
439	>	#endif
440	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
441	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
442	>	endif
443	>	enddo
444	>	if (nGroupTypesRow.eq.0) then
445	>	nGroupTypesRow = nGroupTypesRow + 1
446	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
447	>	groupToGtypeRow(i) = nGroupTypesRow
448	>	else
449	>	GtypeFound = .false.
450	>	do g = 1, nGroupTypesRow
451	>	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
452	>	groupToGtypeRow(i) = g
453	>	GtypeFound = .true.
454	>	endif
455	>	enddo
456	>	if (.not.GtypeFound) then
457	>	nGroupTypesRow = nGroupTypesRow + 1
458	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
459	>	groupToGtypeRow(i) = nGroupTypesRow
460		endif
461	+	endif
462	+	enddo
463	+
464	+	#ifdef IS_MPI
465	+	do j = jstart, jend
466	+	n_in_j = groupStartCol(j+1) - groupStartCol(j)
467	+	groupMaxCutoffCol(j) = 0.0_dp
468	+	do ja = groupStartCol(j), groupStartCol(j+1)-1
469	+	atom1 = groupListCol(ja)
470	+
471	+	me_j = atid_col(atom1)
472	+
473	+	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
474	+	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
475	+	endif
476	+	enddo
477	+
478	+	if (nGroupTypesCol.eq.0) then
479	+	nGroupTypesCol = nGroupTypesCol + 1
480	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
481	+	groupToGtypeCol(j) = nGroupTypesCol
482	+	else
483	+	GtypeFound = .false.
484	+	do g = 1, nGroupTypesCol
485	+	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
486	+	groupToGtypeCol(j) = g
487	+	GtypeFound = .true.
488	+	endif
489	+	enddo
490	+	if (.not.GtypeFound) then
491	+	nGroupTypesCol = nGroupTypesCol + 1
492	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
493	+	groupToGtypeCol(j) = nGroupTypesCol
494	+	endif
495	+	endif
496	+	enddo
497	+
498	+	#else
499	+	! Set pointers to information we just found
500	+	nGroupTypesCol = nGroupTypesRow
501	+	groupToGtypeCol => groupToGtypeRow
502	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
503	+	groupMaxCutoffCol => groupMaxCutoffRow
504	+	#endif
505	+
506	+	!! allocate the gtypeCutoffMap here.
507	+	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
508	+	!! then we do a double loop over all the group TYPES to find the cutoff
509	+	!! map between groups of two types
510	+	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
511	+
512	+	do i = 1, nGroupTypesRow
513	+	do j = 1, nGroupTypesCol
514	+
515	+	select case(cutoffPolicy)
516	+	case(TRADITIONAL_CUTOFF_POLICY)
517	+	thisRcut = tradRcut
518	+	case(MIX_CUTOFF_POLICY)
519	+	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
520	+	case(MAX_CUTOFF_POLICY)
521	+	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
522	+	case default
523	+	call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
524	+	return
525	+	end select
526	+	gtypeCutoffMap(i,j)%rcut = thisRcut
527
528	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
300	<	!             thisRCut = getStickyCutOff(map_i,map_j)
301	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
302	<	endif
303	<
304	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
305	<	!              thisRCut = getStickyPowerCutOff(map_i,map_j)
306	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
307	<	endif
308	<
309	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
310	<	!              thisRCut = getGayberneCutOff(map_i,map_j)
311	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
312	<	endif
313	<
314	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
315	<	!              thisRCut = getGaybrneLJCutOff(map_i,map_j)
316	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
317	<	endif
318	<
319	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
320	<	!              thisRCut = getEAMCutOff(map_i,map_j)
321	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
322	<	endif
323	<
324	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
325	<	!              thisRCut = getShapeCutOff(map_i,map_j)
326	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
327	<	endif
328	<
329	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
330	<	!              thisRCut = getShapeLJCutOff(map_i,map_j)
331	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
332	<	endif
333	<	InteractionMap(map_i, map_j)%rList = actualCutoff
334	<	InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
335	<	end do
336	<	end do
337	<	haveRlist = .true.
338	<	end subroutine createRcuts
528	>	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
529
530	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
531
532	<	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
533	<	!!$  subroutine setRlistDF( this_rlist )
534	<	!!$
344	<	!!$   real(kind=dp) :: this_rlist
345	<	!!$
346	<	!!$    rlist = this_rlist
347	<	!!$    rlistsq = rlist * rlist
348	<	!!$
349	<	!!$    haveRlist = .true.
350	<	!!$
351	<	!!$  end subroutine setRlistDF
532	>	if (.not.haveSkinThickness) then
533	>	skinThickness = 1.0_dp
534	>	endif
535
536	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
537
538	<	subroutine setSimVariables()
355	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
356	<	SIM_uses_LennardJones = SimUsesLennardJones()
357	<	SIM_uses_Electrostatics = SimUsesElectrostatics()
358	<	SIM_uses_Charges = SimUsesCharges()
359	<	SIM_uses_Dipoles = SimUsesDipoles()
360	<	SIM_uses_Sticky = SimUsesSticky()
361	<	SIM_uses_StickyPower = SimUsesStickyPower()
362	<	SIM_uses_GayBerne = SimUsesGayBerne()
363	<	SIM_uses_EAM = SimUsesEAM()
364	<	SIM_uses_Shapes = SimUsesShapes()
365	<	SIM_uses_FLARB = SimUsesFLARB()
366	<	SIM_uses_RF = SimUsesRF()
367	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
368	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
369	<	SIM_uses_PBC = SimUsesPBC()
538	>	! sanity check
539
540	<	haveSIMvariables = .true.
540	>	if (haveDefaultCutoffs) then
541	>	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
542	>	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
543	>	endif
544	>	endif
545	>	enddo
546	>	enddo
547
548	<	return
549	<	end subroutine setSimVariables
548	>	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
549	>	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
550	>	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
551	>	#ifdef IS_MPI
552	>	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
553	>	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
554	>	#endif
555	>	groupMaxCutoffCol => null()
556	>	gtypeMaxCutoffCol => null()
557	>
558	>	haveGtypeCutoffMap = .true.
559	>	end subroutine createGtypeCutoffMap
560
561	+	subroutine setCutoffs(defRcut, defRsw)
562	+
563	+	real(kind=dp),intent(in) :: defRcut, defRsw
564	+	character(len = statusMsgSize) :: errMsg
565	+	integer :: localError
566	+
567	+	defaultRcut = defRcut
568	+	defaultRsw = defRsw
569	+
570	+	defaultDoShift = .false.
571	+	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
572	+
573	+	write(errMsg, *) &
574	+	'cutoffRadius and switchingRadius are set to the same', newline &
575	+	// tab, 'value.  OOPSE will use shifted ', newline &
576	+	// tab, 'potentials instead of switching functions.'
577	+
578	+	call handleInfo("setCutoffs", errMsg)
579	+
580	+	defaultDoShift = .true.
581	+
582	+	endif
583	+
584	+	localError = 0
585	+	call setLJDefaultCutoff( defaultRcut, defaultDoShift )
586	+	call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
587	+	call setCutoffEAM( defaultRcut, localError)
588	+	if (localError /= 0) then
589	+	write(errMsg, *) 'An error has occured in setting the EAM cutoff'
590	+	call handleError("setCutoffs", errMsg)
591	+	end if
592	+	call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
593	+	call setHmatDangerousRcutValue(defaultRcut)
594	+
595	+	haveDefaultCutoffs = .true.
596	+	haveGtypeCutoffMap = .false.
597	+	end subroutine setCutoffs
598	+
599	+	subroutine cWasLame()
600	+
601	+	VisitCutoffsAfterComputing = .true.
602	+	return
603	+
604	+	end subroutine cWasLame
605	+
606	+	subroutine setCutoffPolicy(cutPolicy)
607	+
608	+	integer, intent(in) :: cutPolicy
609	+
610	+	cutoffPolicy = cutPolicy
611	+	haveCutoffPolicy = .true.
612	+	haveGtypeCutoffMap = .false.
613	+
614	+	end subroutine setCutoffPolicy
615	+
616	+	subroutine setElectrostaticMethod( thisESM )
617	+
618	+	integer, intent(in) :: thisESM
619	+
620	+	electrostaticSummationMethod = thisESM
621	+	haveElectrostaticSummationMethod = .true.
622	+
623	+	end subroutine setElectrostaticMethod
624	+
625	+	subroutine setSkinThickness( thisSkin )
626	+
627	+	real(kind=dp), intent(in) :: thisSkin
628	+
629	+	skinThickness = thisSkin
630	+	haveSkinThickness = .true.
631	+	haveGtypeCutoffMap = .false.
632	+
633	+	end subroutine setSkinThickness
634	+
635	+	subroutine setSimVariables()
636	+	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
637	+	SIM_uses_EAM = SimUsesEAM()
638	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
639	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
640	+	SIM_uses_PBC = SimUsesPBC()
641	+	SIM_uses_SC = SimUsesSC()
642	+
643	+	haveSIMvariables = .true.
644	+
645	+	return
646	+	end subroutine setSimVariables
647	+
648		subroutine doReadyCheck(error)
649		integer, intent(out) :: error
650
#	Line 380 | Line 652 | contains
652
653		error = 0
654
655	<	if (.not. haveInteractionMap) then
655	>	if (.not. haveInteractionHash) then
656	>	call createInteractionHash()
657	>	endif
658
659	<	myStatus = 0
659	>	if (.not. haveGtypeCutoffMap) then
660	>	call createGtypeCutoffMap()
661	>	endif
662
387	–	call createInteractionMap(myStatus)
663
664	<	if (myStatus .ne. 0) then
665	<	write(default_error, *) 'createInteractionMap failed in doForces!'
666	<	error = -1
392	<	return
393	<	endif
664	>	if (VisitCutoffsAfterComputing) then
665	>	call set_switch(GROUP_SWITCH, largestRcut, largestRcut)
666	>	call setHmatDangerousRcutValue(largestRcut)
667		endif
668
669	+
670		if (.not. haveSIMvariables) then
671		call setSimVariables()
672		endif
673
674	<	if (.not. haveRlist) then
675	<	write(default_error, *) 'rList has not been set in doForces!'
676	<	error = -1
677	<	return
678	<	endif
674	>	!  if (.not. haveRlist) then
675	>	!     write(default_error, *) 'rList has not been set in doForces!'
676	>	!     error = -1
677	>	!     return
678	>	!  endif
679
680		if (.not. haveNeighborList) then
681		write(default_error, *) 'neighbor list has not been initialized in doForces!'
#	Line 426 | Line 700 | contains
700		end subroutine doReadyCheck
701
702
703	<	subroutine init_FF(use_RF_c, thisStat)
703	>	subroutine init_FF(thisStat)
704
431	–	logical, intent(in) :: use_RF_c
432	–
705		integer, intent(out) :: thisStat
706		integer :: my_status, nMatches
707		integer, pointer :: MatchList(:) => null()
436	–	real(kind=dp) :: rcut, rrf, rt, dielect
708
709		!! assume things are copacetic, unless they aren't
710		thisStat = 0
711
441	–	!! Fortran's version of a cast:
442	–	FF_uses_RF = use_RF_c
443	–
712		!! init_FF is called *after* all of the atom types have been
713		!! defined in atype_module using the new_atype subroutine.
714		!!
#	Line 448 | Line 716 | contains
716		!! interactions are used by the force field.
717
718		FF_uses_DirectionalAtoms = .false.
451	–	FF_uses_LennardJones = .false.
452	–	FF_uses_Electrostatics = .false.
453	–	FF_uses_Charges = .false.
719		FF_uses_Dipoles = .false.
455	–	FF_uses_Sticky = .false.
456	–	FF_uses_StickyPower = .false.
720		FF_uses_GayBerne = .false.
721		FF_uses_EAM = .false.
722	<	FF_uses_Shapes = .false.
460	<	FF_uses_FLARB = .false.
722	>	FF_uses_SC = .false.
723
724		call getMatchingElementList(atypes, "is_Directional", .true., &
725		nMatches, MatchList)
726		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
727
466	–	call getMatchingElementList(atypes, "is_LennardJones", .true., &
467	–	nMatches, MatchList)
468	–	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
469	–
470	–	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
471	–	nMatches, MatchList)
472	–	if (nMatches .gt. 0) then
473	–	FF_uses_Electrostatics = .true.
474	–	endif
475	–
476	–	call getMatchingElementList(atypes, "is_Charge", .true., &
477	–	nMatches, MatchList)
478	–	if (nMatches .gt. 0) then
479	–	FF_uses_Charges = .true.
480	–	FF_uses_Electrostatics = .true.
481	–	endif
482	–
728		call getMatchingElementList(atypes, "is_Dipole", .true., &
729		nMatches, MatchList)
730	<	if (nMatches .gt. 0) then
486	<	FF_uses_Dipoles = .true.
487	<	FF_uses_Electrostatics = .true.
488	<	FF_uses_DirectionalAtoms = .true.
489	<	endif
490	<
491	<	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
492	<	nMatches, MatchList)
493	<	if (nMatches .gt. 0) then
494	<	FF_uses_Quadrupoles = .true.
495	<	FF_uses_Electrostatics = .true.
496	<	FF_uses_DirectionalAtoms = .true.
497	<	endif
498	<
499	<	call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
500	<	MatchList)
501	<	if (nMatches .gt. 0) then
502	<	FF_uses_Sticky = .true.
503	<	FF_uses_DirectionalAtoms = .true.
504	<	endif
505	<
506	<	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
507	<	MatchList)
508	<	if (nMatches .gt. 0) then
509	<	FF_uses_StickyPower = .true.
510	<	FF_uses_DirectionalAtoms = .true.
511	<	endif
730	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
731
732		call getMatchingElementList(atypes, "is_GayBerne", .true., &
733		nMatches, MatchList)
734	<	if (nMatches .gt. 0) then
516	<	FF_uses_GayBerne = .true.
517	<	FF_uses_DirectionalAtoms = .true.
518	<	endif
734	>	if (nMatches .gt. 0) FF_uses_GayBerne = .true.
735
736		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
737		if (nMatches .gt. 0) FF_uses_EAM = .true.
738
739	<	call getMatchingElementList(atypes, "is_Shape", .true., &
740	<	nMatches, MatchList)
525	<	if (nMatches .gt. 0) then
526	<	FF_uses_Shapes = .true.
527	<	FF_uses_DirectionalAtoms = .true.
528	<	endif
739	>	call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
740	>	if (nMatches .gt. 0) FF_uses_SC = .true.
741
530	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
531	–	nMatches, MatchList)
532	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
742
534	–	!! Assume sanity (for the sake of argument)
743		haveSaneForceField = .true.
536	–
537	–	!! check to make sure the FF_uses_RF setting makes sense
538	–
539	–	if (FF_uses_dipoles) then
540	–	if (FF_uses_RF) then
541	–	dielect = getDielect()
542	–	call initialize_rf(dielect)
543	–	endif
544	–	else
545	–	if (FF_uses_RF) then
546	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
547	–	thisStat = -1
548	–	haveSaneForceField = .false.
549	–	return
550	–	endif
551	–	endif
744
553	–	!sticky module does not contain check_sticky_FF anymore
554	–	!if (FF_uses_sticky) then
555	–	!   call check_sticky_FF(my_status)
556	–	!   if (my_status /= 0) then
557	–	!      thisStat = -1
558	–	!      haveSaneForceField = .false.
559	–	!      return
560	–	!   end if
561	–	!endif
562	–
745		if (FF_uses_EAM) then
746		call init_EAM_FF(my_status)
747		if (my_status /= 0) then
#	Line 570 | Line 752 | contains
752		end if
753		endif
754
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	–
582	–	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
583	–	endif
584	–
755		if (.not. haveNeighborList) then
756		!! Create neighbor lists
757		call expandNeighborList(nLocal, my_status)
#	Line 615 | Line 785 | contains
785
786		!! Stress Tensor
787		real( kind = dp), dimension(9) :: tau
788	<	real ( kind = dp ) :: pot
788	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
789		logical ( kind = 2) :: do_pot_c, do_stress_c
790		logical :: do_pot
791		logical :: do_stress
792		logical :: in_switching_region
793		#ifdef IS_MPI
794	<	real( kind = DP ) :: pot_local
794	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
795		integer :: nAtomsInRow
796		integer :: nAtomsInCol
797		integer :: nprocs
#	Line 636 | Line 806 | contains
806		integer :: nlist
807		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
808		real( kind = DP ) :: sw, dswdr, swderiv, mf
809	+	real( kind = DP ) :: rVal
810		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
811		real(kind=dp) :: rfpot, mu_i, virial
812	+	real(kind=dp):: rCut
813		integer :: me_i, me_j, n_in_i, n_in_j
814		logical :: is_dp_i
815		integer :: neighborListSize
#	Line 645 | Line 817 | contains
817		integer :: localError
818		integer :: propPack_i, propPack_j
819		integer :: loopStart, loopEnd, loop
820	<	integer :: iMap
821	<	real(kind=dp) :: listSkin = 1.0
820	>	integer :: iHash
821	>	integer :: i1
822	>
823
824		!! initialize local variables
825
#	Line 710 | Line 883 | contains
883		! (but only on the first time through):
884		if (loop .eq. loopStart) then
885		#ifdef IS_MPI
886	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
886	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
887		update_nlist)
888		#else
889	<	call checkNeighborList(nGroups, q_group, listSkin, &
889	>	call checkNeighborList(nGroups, q_group, skinThickness, &
890		update_nlist)
891		#endif
892		endif
#	Line 743 | Line 916 | contains
916
917		if (update_nlist) then
918		#ifdef IS_MPI
746	–	me_i = atid_row(i)
919		jstart = 1
920		jend = nGroupsInCol
921		#else
750	–	me_i = atid(i)
922		jstart = i+1
923		jend = nGroups
924		#endif
#	Line 773 | Line 944 | contains
944		me_j = atid(j)
945		call get_interatomic_vector(q_group(:,i), &
946		q_group(:,j), d_grp, rgrpsq)
947	<	#endif
947	>	#endif
948
949	<	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
949	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
950		if (update_nlist) then
951		nlist = nlist + 1
952
#	Line 795 | Line 966 | contains
966
967		list(nlist) = j
968		endif
969	+
970
971	<	if (loop .eq. PAIR_LOOP) then
972	<	vij = 0.0d0
801	<	fij(1:3) = 0.0d0
802	<	endif
971	>
972	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
973
974	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
975	<	in_switching_region)
976	<
977	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
978	<
979	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
980	<
981	<	atom1 = groupListRow(ia)
982	<
983	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
984	<
985	<	atom2 = groupListCol(jb)
986	<
987	<	if (skipThisPair(atom1, atom2)) cycle inner
988	<
989	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
990	<	d_atm(1:3) = d_grp(1:3)
991	<	ratmsq = rgrpsq
992	<	else
993	<	#ifdef IS_MPI
994	<	call get_interatomic_vector(q_Row(:,atom1), &
995	<	q_Col(:,atom2), d_atm, ratmsq)
974	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
975	>	if (loop .eq. PAIR_LOOP) then
976	>	vij = 0.0d0
977	>	fij(1:3) = 0.0d0
978	>	endif
979	>
980	>	call get_switch(rgrpsq, sw, dswdr, rgrp, &
981	>	group_switch, in_switching_region)
982	>
983	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
984	>
985	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
986	>
987	>	atom1 = groupListRow(ia)
988	>
989	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
990	>
991	>	atom2 = groupListCol(jb)
992	>
993	>	if (skipThisPair(atom1, atom2))  cycle inner
994	>
995	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
996	>	d_atm(1:3) = d_grp(1:3)
997	>	ratmsq = rgrpsq
998	>	else
999	>	#ifdef IS_MPI
1000	>	call get_interatomic_vector(q_Row(:,atom1), &
1001	>	q_Col(:,atom2), d_atm, ratmsq)
1002		#else
1003	<	call get_interatomic_vector(q(:,atom1), &
1004	<	q(:,atom2), d_atm, ratmsq)
1003	>	call get_interatomic_vector(q(:,atom1), &
1004	>	q(:,atom2), d_atm, ratmsq)
1005		#endif
1006	<	endif
1007	<
1008	<	if (loop .eq. PREPAIR_LOOP) then
1006	>	endif
1007	>
1008	>	if (loop .eq. PREPAIR_LOOP) then
1009		#ifdef IS_MPI
1010	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1011	<	rgrpsq, d_grp, do_pot, do_stress, &
1012	<	eFrame, A, f, t, pot_local)
1010	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1011	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1012	>	eFrame, A, f, t, pot_local)
1013		#else
1014	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1015	<	rgrpsq, d_grp, do_pot, do_stress, &
1016	<	eFrame, A, f, t, pot)
1014	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1015	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1016	>	eFrame, A, f, t, pot)
1017		#endif
1018	<	else
1018	>	else
1019		#ifdef IS_MPI
1020	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1021	<	do_pot, &
1022	<	eFrame, A, f, t, pot_local, vpair, fpair)
1020	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1021	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1022	>	fpair, d_grp, rgrp, rCut)
1023		#else
1024	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1025	<	do_pot,  &
1026	<	eFrame, A, f, t, pot, vpair, fpair)
1024	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1025	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1026	>	d_grp, rgrp, rCut)
1027		#endif
1028	+	vij = vij + vpair
1029	+	fij(1:3) = fij(1:3) + fpair(1:3)
1030	+	endif
1031	+	enddo inner
1032	+	enddo
1033
1034	<	vij = vij + vpair
1035	<	fij(1:3) = fij(1:3) + fpair(1:3)
1036	<	endif
1037	<	enddo inner
1038	<	enddo
1039	<
1040	<	if (loop .eq. PAIR_LOOP) then
1041	<	if (in_switching_region) then
1042	<	swderiv = vij*dswdr/rgrp
1043	<	fij(1) = fij(1) + swderiv*d_grp(1)
863	<	fij(2) = fij(2) + swderiv*d_grp(2)
864	<	fij(3) = fij(3) + swderiv*d_grp(3)
865	<
866	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
867	<	atom1=groupListRow(ia)
868	<	mf = mfactRow(atom1)
1034	>	if (loop .eq. PAIR_LOOP) then
1035	>	if (in_switching_region) then
1036	>	swderiv = vij*dswdr/rgrp
1037	>	fij(1) = fij(1) + swderiv*d_grp(1)
1038	>	fij(2) = fij(2) + swderiv*d_grp(2)
1039	>	fij(3) = fij(3) + swderiv*d_grp(3)
1040	>
1041	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1042	>	atom1=groupListRow(ia)
1043	>	mf = mfactRow(atom1)
1044		#ifdef IS_MPI
1045	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1046	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1047	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1045	>	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1046	>	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1047	>	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1048		#else
1049	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1050	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1051	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1049	>	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1050	>	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1051	>	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1052		#endif
1053	<	enddo
1054	<
1055	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1056	<	atom2=groupListCol(jb)
1057	<	mf = mfactCol(atom2)
1053	>	enddo
1054	>
1055	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1056	>	atom2=groupListCol(jb)
1057	>	mf = mfactCol(atom2)
1058		#ifdef IS_MPI
1059	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1060	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1061	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1059	>	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1060	>	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1061	>	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1062		#else
1063	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1064	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1065	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1063	>	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1064	>	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1065	>	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1066		#endif
1067	<	enddo
1068	<	endif
1067	>	enddo
1068	>	endif
1069
1070	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1070	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1071	>	endif
1072		endif
1073	<	end if
1073	>	endif
1074		enddo
1075	+
1076		enddo outer
1077
1078		if (update_nlist) then
#	Line 955 | Line 1132 | contains
1132
1133		if (do_pot) then
1134		! scatter/gather pot_row into the members of my column
1135	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1136	<
1135	>	do i = 1,LR_POT_TYPES
1136	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1137	>	end do
1138		! scatter/gather pot_local into all other procs
1139		! add resultant to get total pot
1140		do i = 1, nlocal
1141	<	pot_local = pot_local + pot_Temp(i)
1141	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1142	>	+ pot_Temp(1:LR_POT_TYPES,i)
1143		enddo
1144
1145		pot_Temp = 0.0_DP
1146	<
1147	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1146	>	do i = 1,LR_POT_TYPES
1147	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1148	>	end do
1149		do i = 1, nlocal
1150	<	pot_local = pot_local + pot_Temp(i)
1150	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1151	>	+ pot_Temp(1:LR_POT_TYPES,i)
1152		enddo
1153
1154		endif
1155		#endif
1156
1157	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1157	>	if (SIM_requires_postpair_calc) then
1158	>	do i = 1, nlocal
1159	>
1160	>	! we loop only over the local atoms, so we don't need row and column
1161	>	! lookups for the types
1162	>
1163	>	me_i = atid(i)
1164	>
1165	>	! is the atom electrostatic?  See if it would have an
1166	>	! electrostatic interaction with itself
1167	>	iHash = InteractionHash(me_i,me_i)
1168
1169	<	if (FF_uses_RF .and. SIM_uses_RF) then
979	<
1169	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1170		#ifdef IS_MPI
1171	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1172	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
983	<	do i = 1,nlocal
984	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
985	<	end do
986	<	#endif
987	<
988	<	do i = 1, nLocal
989	<
990	<	rfpot = 0.0_DP
991	<	#ifdef IS_MPI
992	<	me_i = atid_row(i)
1171	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1172	>	t, do_pot)
1173		#else
1174	<	me_i = atid(i)
1174	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1175	>	t, do_pot)
1176		#endif
1177	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1177	>	endif
1178	>
1179	>
1180	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1181
1182	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1183	<
1184	<	mu_i = getDipoleMoment(me_i)
1185	<
1186	<	!! The reaction field needs to include a self contribution
1187	<	!! to the field:
1188	<	call accumulate_self_rf(i, mu_i, eFrame)
1189	<	!! Get the reaction field contribution to the
1190	<	!! potential and torques:
1191	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1182	>	! loop over the excludes to accumulate RF stuff we've
1183	>	! left out of the normal pair loop
1184	>
1185	>	do i1 = 1, nSkipsForAtom(i)
1186	>	j = skipsForAtom(i, i1)
1187	>
1188	>	! prevent overcounting of the skips
1189	>	if (i.lt.j) then
1190	>	call get_interatomic_vector(q(:,i), &
1191	>	q(:,j), d_atm, ratmsq)
1192	>	rVal = dsqrt(ratmsq)
1193	>	call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1194	>	in_switching_region)
1195		#ifdef IS_MPI
1196	<	pot_local = pot_local + rfpot
1196	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1197	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1198		#else
1199	<	pot = pot + rfpot
1200	<
1199	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1200	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1201		#endif
1202	<	endif
1203	<	enddo
1204	<	endif
1202	>	endif
1203	>	enddo
1204	>	endif
1205	>	enddo
1206		endif
1207	<
1019	<
1207	>
1208		#ifdef IS_MPI
1209	<
1209	>
1210		if (do_pot) then
1211	<	pot = pot + pot_local
1212	<	!! we assume the c code will do the allreduce to get the total potential
1025	<	!! we could do it right here if we needed to...
1211	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1212	>	mpi_comm_world,mpi_err)
1213		endif
1214	<
1214	>
1215		if (do_stress) then
1216		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1217		mpi_comm_world,mpi_err)
1218		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1219		mpi_comm_world,mpi_err)
1220		endif
1221	<
1221	>
1222		#else
1223	<
1223	>
1224		if (do_stress) then
1225		tau = tau_Temp
1226		virial = virial_Temp
1227		endif
1228	<
1228	>
1229		#endif
1230	<
1230	>
1231		end subroutine do_force_loop
1232
1233		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1234	<	eFrame, A, f, t, pot, vpair, fpair)
1234	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1235
1236	<	real( kind = dp ) :: pot, vpair, sw
1236	>	real( kind = dp ) :: vpair, sw
1237	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1238		real( kind = dp ), dimension(3) :: fpair
1239		real( kind = dp ), dimension(nLocal)   :: mfact
1240		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1057 | Line 1245 | contains
1245		logical, intent(inout) :: do_pot
1246		integer, intent(in) :: i, j
1247		real ( kind = dp ), intent(inout) :: rijsq
1248	<	real ( kind = dp )                :: r
1248	>	real ( kind = dp ), intent(inout) :: r_grp
1249		real ( kind = dp ), intent(inout) :: d(3)
1250	<	real ( kind = dp ) :: ebalance
1250	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1251	>	real ( kind = dp ), intent(inout) :: rCut
1252	>	real ( kind = dp ) :: r
1253		integer :: me_i, me_j
1254
1255	<	integer :: iMap
1255	>	integer :: iHash
1256
1257		r = sqrt(rijsq)
1258		vpair = 0.0d0
#	Line 1076 | Line 1266 | contains
1266		me_j = atid(j)
1267		#endif
1268
1269	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1270	<
1271	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1272	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1273	<	endif
1084	<
1085	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1086	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1087	<	pot, eFrame, f, t, do_pot)
1088	<
1089	<	if (FF_uses_RF .and. SIM_uses_RF) then
1090	<
1091	<	! CHECK ME (RF needs to know about all electrostatic types)
1092	<	call accumulate_rf(i, j, r, eFrame, sw)
1093	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1094	<	endif
1095	<
1269	>	iHash = InteractionHash(me_i, me_j)
1270	>
1271	>	if ( iand(iHash, LJ_PAIR).ne.0 ) then
1272	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1273	>	pot(VDW_POT), f, do_pot)
1274		endif
1275	<
1276	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1275	>
1276	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1277	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1278	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1279	>	endif
1280	>
1281	>	if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1282		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1283	<	pot, A, f, t, do_pot)
1283	>	pot(HB_POT), A, f, t, do_pot)
1284		endif
1285	<
1286	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1285	>
1286	>	if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1287		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1288	<	pot, A, f, t, do_pot)
1288	>	pot(HB_POT), A, f, t, do_pot)
1289		endif
1290	<
1291	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1290	>
1291	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1292		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1293	<	pot, A, f, t, do_pot)
1293	>	pot(VDW_POT), A, f, t, do_pot)
1294		endif
1295
1296	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1297	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1298	<	!           pot, A, f, t, do_pot)
1296	>	if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1297	>	call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1298	>	pot(VDW_POT), A, f, t, do_pot)
1299		endif
1300	<
1301	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1302	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1303	<	do_pot)
1300	>
1301	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1302	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1303	>	pot(METALLIC_POT), f, do_pot)
1304		endif
1305	<
1306	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1305	>
1306	>	if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1307		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1308	<	pot, A, f, t, do_pot)
1308	>	pot(VDW_POT), A, f, t, do_pot)
1309		endif
1310	<
1311	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1310	>
1311	>	if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1312		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1313	<	pot, A, f, t, do_pot)
1313	>	pot(VDW_POT), A, f, t, do_pot)
1314		endif
1315	+
1316	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1317	+	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1318	+	pot(METALLIC_POT), f, do_pot)
1319	+	endif
1320	+
1321
1322	+
1323		end subroutine do_pair
1324
1325	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1325	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1326		do_pot, do_stress, eFrame, A, f, t, pot)
1327
1328	<	real( kind = dp ) :: pot, sw
1328	>	real( kind = dp ) :: sw
1329	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1330		real( kind = dp ), dimension(9,nLocal) :: eFrame
1331		real (kind=dp), dimension(9,nLocal) :: A
1332		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1143 | Line 1334 | contains
1334
1335		logical, intent(inout) :: do_pot, do_stress
1336		integer, intent(in) :: i, j
1337	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1337	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1338		real ( kind = dp )                :: r, rc
1339		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1340
1341	<	integer :: me_i, me_j, iMap
1341	>	integer :: me_i, me_j, iHash
1342
1343	+	r = sqrt(rijsq)
1344	+
1345		#ifdef IS_MPI
1346		me_i = atid_row(i)
1347		me_j = atid_col(j)
#	Line 1157 | Line 1350 | contains
1350		me_j = atid(j)
1351		#endif
1352
1353	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1353	>	iHash = InteractionHash(me_i, me_j)
1354
1355	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1356	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1355	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1356	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1357		endif
1358	+
1359	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1360	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1361	+	endif
1362
1363		end subroutine do_prepair
1364
1365
1366		subroutine do_preforce(nlocal,pot)
1367		integer :: nlocal
1368	<	real( kind = dp ) :: pot
1368	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1369
1370		if (FF_uses_EAM .and. SIM_uses_EAM) then
1371	<	call calc_EAM_preforce_Frho(nlocal,pot)
1371	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1372		endif
1373	+	if (FF_uses_SC .and. SIM_uses_SC) then
1374	+	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1375	+	endif
1376
1377
1378		end subroutine do_preforce
#	Line 1257 | Line 1457 | contains
1457		pot_Col = 0.0_dp
1458		pot_Temp = 0.0_dp
1459
1260	–	rf_Row = 0.0_dp
1261	–	rf_Col = 0.0_dp
1262	–	rf_Temp = 0.0_dp
1263	–
1460		#endif
1461
1462		if (FF_uses_EAM .and. SIM_uses_EAM) then
1463		call clean_EAM()
1464		endif
1465
1270	–	rf = 0.0_dp
1466		tau_Temp = 0.0_dp
1467		virial_Temp = 0.0_dp
1468		end subroutine zero_work_arrays
#	Line 1356 | Line 1551 | contains
1551
1552		function FF_UsesDirectionalAtoms() result(doesit)
1553		logical :: doesit
1554	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1360	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1361	<	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1554	>	doesit = FF_uses_DirectionalAtoms
1555		end function FF_UsesDirectionalAtoms
1556
1557		function FF_RequiresPrepairCalc() result(doesit)
1558		logical :: doesit
1559	<	doesit = FF_uses_EAM
1559	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1560	>	.or. FF_uses_MEAM
1561		end function FF_RequiresPrepairCalc
1562
1369	–	function FF_RequiresPostpairCalc() result(doesit)
1370	–	logical :: doesit
1371	–	doesit = FF_uses_RF
1372	–	end function FF_RequiresPostpairCalc
1373	–
1563		#ifdef PROFILE
1564		function getforcetime() result(totalforcetime)
1565		real(kind=dp) :: totalforcetime

Diff Legend

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