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 2267 by tim, Fri Jul 29 17:34:06 2005 UTC vs.
Revision 2715 by chrisfen, Sun Apr 16 02:51:16 2006 UTC

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

Diff Legend

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