[ViewVC] Diff of: group/trunk/OOPSE-4/src/UseTheForce/doForces.F90

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2261 by gezelter, Tue Jun 28 13:58:45 2005 UTC vs.
Revision 2411 by chrisfen, Wed Nov 2 21:01:21 2005 UTC

#	Line 45 \| Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.22 2005-06-28 13:58:45 gezelter Exp $, $Date: 2005-06-28 13:58:45 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $
48	>	!! @version $Id: doForces.F90,v 1.67 2005-11-02 21:01:18 chrisfen Exp $, $Date: 2005-11-02 21:01:18 $, $Name: not supported by cvs2svn $, $Revision: 1.67 $
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
#	Line 73 \| Line 72 \| module doForces
72
73		#define __FORTRAN90
74		#include "UseTheForce/fSwitchingFunction.h"
75	+	#include "UseTheForce/fCutoffPolicy.h"
76		#include "UseTheForce/DarkSide/fInteractionMap.h"
77	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78
79	+
80		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81		INTEGER, PARAMETER:: PAIR_LOOP = 2
82
81	–	logical, save :: haveRlist = .false.
83		logical, save :: haveNeighborList = .false.
84		logical, save :: haveSIMvariables = .false.
85		logical, save :: haveSaneForceField = .false.
86	+	logical, save :: haveInteractionHash = .false.
87	+	logical, save :: haveGtypeCutoffMap = .false.
88	+	logical, save :: haveDefaultCutoffs = .false.
89	+	logical, save :: haveRlist = .false.
90
91		logical, save :: FF_uses_DirectionalAtoms
87	–	logical, save :: FF_uses_LennardJones
88	–	logical, save :: FF_uses_Electrostatics
89	–	logical, save :: FF_uses_Charges
92		logical, save :: FF_uses_Dipoles
91	–	logical, save :: FF_uses_Quadrupoles
92	–	logical, save :: FF_uses_Sticky
93	–	logical, save :: FF_uses_StickyPower
93		logical, save :: FF_uses_GayBerne
94		logical, save :: FF_uses_EAM
96	–	logical, save :: FF_uses_Shapes
97	–	logical, save :: FF_uses_FLARB
98	–	logical, save :: FF_uses_RF
95
96		logical, save :: SIM_uses_DirectionalAtoms
101	–	logical, save :: SIM_uses_LennardJones
102	–	logical, save :: SIM_uses_Electrostatics
103	–	logical, save :: SIM_uses_Charges
104	–	logical, save :: SIM_uses_Dipoles
105	–	logical, save :: SIM_uses_Quadrupoles
106	–	logical, save :: SIM_uses_Sticky
107	–	logical, save :: SIM_uses_StickyPower
108	–	logical, save :: SIM_uses_GayBerne
97		logical, save :: SIM_uses_EAM
110	–	logical, save :: SIM_uses_Shapes
111	–	logical, save :: SIM_uses_FLARB
112	–	logical, save :: SIM_uses_RF
98		logical, save :: SIM_requires_postpair_calc
99		logical, save :: SIM_requires_prepair_calc
100		logical, save :: SIM_uses_PBC
116	–	logical, save :: SIM_uses_molecular_cutoffs
101
102	<	!!!GO AWAY---------
119	<	!!!!!real(kind=dp), save :: rlist, rlistsq
102	>	integer, save :: electrostaticSummationMethod
103
104		public :: init_FF
105	+	public :: setDefaultCutoffs
106		public :: do_force_loop
107	<	! public :: setRlistDF
107	>	public :: createInteractionHash
108	>	public :: createGtypeCutoffMap
109	>	public :: getStickyCut
110	>	public :: getStickyPowerCut
111	>	public :: getGayBerneCut
112	>	public :: getEAMCut
113	>	public :: getShapeCut
114
115		#ifdef PROFILE
116		public :: getforcetime
#	Line 128 \| Line 118 \| module doForces
118		real :: forceTimeInitial, forceTimeFinal
119		integer :: nLoops
120		#endif
131	–
132	–	type, public :: Interaction
133	–	integer :: InteractionHash
134	–	real(kind=dp) :: rCut
135	–	end type Interaction
121
122	<	type(Interaction), dimension(:,:),allocatable :: InteractionMap
123	<
124	<	!public :: addInteraction
125	<	!public :: setInteractionHash
126	<	!public :: getInteractionHash
127	<	public :: createInteractionMap
122	>	!! Variables for cutoff mapping and interaction mapping
123	>	! Bit hash to determine pair-pair interactions.
124	>	integer, dimension(:,:), allocatable :: InteractionHash
125	>	real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
126	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
127	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
128	>
129	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
130	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
131	>
132	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
133	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
134	>	type ::gtypeCutoffs
135	>	real(kind=dp) :: rcut
136	>	real(kind=dp) :: rcutsq
137	>	real(kind=dp) :: rlistsq
138	>	end type gtypeCutoffs
139	>	type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
140	>
141	>	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
142	>	real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
143	>	real(kind=dp),save :: listSkin
144
145		contains
146
147	<
147	<	subroutine createInteractionMap(status)
147	>	subroutine createInteractionHash(status)
148		integer :: nAtypes
149	<	integer :: status
149	>	integer, intent(out) :: status
150		integer :: i
151		integer :: j
152	<	integer :: ihash
153	<	real(kind=dp) :: myRcut
154	<	! Test Types
152	>	integer :: iHash
153	>	!! Test Types
154		logical :: i_is_LJ
155		logical :: i_is_Elect
156		logical :: i_is_Sticky
#	Line 166 \| Line 165 \| contains
165		logical :: j_is_GB
166		logical :: j_is_EAM
167		logical :: j_is_Shape
168	<
169	<
168	>	real(kind=dp) :: myRcut
169	>
170	>	status = 0
171	>
172		if (.not. associated(atypes)) then
173	<	call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
173	>	call handleError("atype", "atypes was not present before call of createInteractionHash!")
174		status = -1
175		return
176		endif
#	Line 181 \| Line 182 \| contains
182		return
183		end if
184
185	<	if (.not. allocated(InteractionMap)) then
186	<	allocate(InteractionMap(nAtypes,nAtypes))
185	>	if (.not. allocated(InteractionHash)) then
186	>	allocate(InteractionHash(nAtypes,nAtypes))
187	>	else
188	>	deallocate(InteractionHash)
189	>	allocate(InteractionHash(nAtypes,nAtypes))
190		endif
191	+
192	+	if (.not. allocated(atypeMaxCutoff)) then
193	+	allocate(atypeMaxCutoff(nAtypes))
194	+	else
195	+	deallocate(atypeMaxCutoff)
196	+	allocate(atypeMaxCutoff(nAtypes))
197	+	endif
198
199		do i = 1, nAtypes
200		call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
#	Line 236 \| Line 247 \| contains
247		if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
248
249
250	<	InteractionMap(i,j)%InteractionHash = iHash
251	<	InteractionMap(j,i)%InteractionHash = iHash
250	>	InteractionHash(i,j) = iHash
251	>	InteractionHash(j,i) = iHash
252
253		end do
254
255		end do
245	–	end subroutine createInteractionMap
256
257	+	haveInteractionHash = .true.
258	+	end subroutine createInteractionHash
259
260	+	subroutine createGtypeCutoffMap(stat)
261
262	<	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
263	<	!!$ subroutine setRlistDF( this_rlist )
264	<	!!$
265	<	!!$ real(kind=dp) :: this_rlist
266	<	!!$
267	<	!!$ rlist = this_rlist
268	<	!!$ rlistsq = rlist * rlist
269	<	!!$
270	<	!!$ haveRlist = .true.
258	<	!!$
259	<	!!$ end subroutine setRlistDF
262	>	integer, intent(out), optional :: stat
263	>	logical :: i_is_LJ
264	>	logical :: i_is_Elect
265	>	logical :: i_is_Sticky
266	>	logical :: i_is_StickyP
267	>	logical :: i_is_GB
268	>	logical :: i_is_EAM
269	>	logical :: i_is_Shape
270	>	logical :: GtypeFound
271
272	+	integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend
273	+	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
274	+	integer :: nGroupsInRow
275	+	integer :: nGroupsInCol
276	+	integer :: nGroupTypesRow,nGroupTypesCol
277	+	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
278	+	real(kind=dp) :: biggestAtypeCutoff
279
280	+	stat = 0
281	+	if (.not. haveInteractionHash) then
282	+	call createInteractionHash(myStatus)
283	+	if (myStatus .ne. 0) then
284	+	write(default_error, *) 'createInteractionHash failed in doForces!'
285	+	stat = -1
286	+	return
287	+	endif
288	+	endif
289	+	#ifdef IS_MPI
290	+	nGroupsInRow = getNgroupsInRow(plan_group_row)
291	+	nGroupsInCol = getNgroupsInCol(plan_group_col)
292	+	#endif
293	+	nAtypes = getSize(atypes)
294	+	! Set all of the initial cutoffs to zero.
295	+	atypeMaxCutoff = 0.0_dp
296	+	do i = 1, nAtypes
297	+	if (SimHasAtype(i)) then
298	+	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
299	+	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
300	+	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
301	+	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
302	+	call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
303	+	call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
304	+	call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
305	+
306	+
307	+	if (haveDefaultCutoffs) then
308	+	atypeMaxCutoff(i) = defaultRcut
309	+	else
310	+	if (i_is_LJ) then
311	+	thisRcut = getSigma(i) * 2.5_dp
312	+	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313	+	endif
314	+	if (i_is_Elect) then
315	+	thisRcut = defaultRcut
316	+	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317	+	endif
318	+	if (i_is_Sticky) then
319	+	thisRcut = getStickyCut(i)
320	+	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321	+	endif
322	+	if (i_is_StickyP) then
323	+	thisRcut = getStickyPowerCut(i)
324	+	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325	+	endif
326	+	if (i_is_GB) then
327	+	thisRcut = getGayBerneCut(i)
328	+	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329	+	endif
330	+	if (i_is_EAM) then
331	+	thisRcut = getEAMCut(i)
332	+	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333	+	endif
334	+	if (i_is_Shape) then
335	+	thisRcut = getShapeCut(i)
336	+	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337	+	endif
338	+	endif
339	+
340	+
341	+	if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
342	+	biggestAtypeCutoff = atypeMaxCutoff(i)
343	+	endif
344	+
345	+	endif
346	+	enddo
347	+
348	+
349	+
350	+	istart = 1
351	+	jstart = 1
352	+	#ifdef IS_MPI
353	+	iend = nGroupsInRow
354	+	jend = nGroupsInCol
355	+	#else
356	+	iend = nGroups
357	+	jend = nGroups
358	+	#endif
359	+
360	+	!! allocate the groupToGtype and gtypeMaxCutoff here.
361	+	if(.not.allocated(groupToGtypeRow)) then
362	+	! allocate(groupToGtype(iend))
363	+	allocate(groupToGtypeRow(iend))
364	+	else
365	+	deallocate(groupToGtypeRow)
366	+	allocate(groupToGtypeRow(iend))
367	+	endif
368	+	if(.not.allocated(groupMaxCutoffRow)) then
369	+	allocate(groupMaxCutoffRow(iend))
370	+	else
371	+	deallocate(groupMaxCutoffRow)
372	+	allocate(groupMaxCutoffRow(iend))
373	+	end if
374	+
375	+	if(.not.allocated(gtypeMaxCutoffRow)) then
376	+	allocate(gtypeMaxCutoffRow(iend))
377	+	else
378	+	deallocate(gtypeMaxCutoffRow)
379	+	allocate(gtypeMaxCutoffRow(iend))
380	+	endif
381	+
382	+
383	+	#ifdef IS_MPI
384	+	! We only allocate new storage if we are in MPI because Ncol /= Nrow
385	+	if(.not.associated(groupToGtypeCol)) then
386	+	allocate(groupToGtypeCol(jend))
387	+	else
388	+	deallocate(groupToGtypeCol)
389	+	allocate(groupToGtypeCol(jend))
390	+	end if
391	+
392	+	if(.not.associated(groupToGtypeCol)) then
393	+	allocate(groupToGtypeCol(jend))
394	+	else
395	+	deallocate(groupToGtypeCol)
396	+	allocate(groupToGtypeCol(jend))
397	+	end if
398	+	if(.not.associated(gtypeMaxCutoffCol)) then
399	+	allocate(gtypeMaxCutoffCol(jend))
400	+	else
401	+	deallocate(gtypeMaxCutoffCol)
402	+	allocate(gtypeMaxCutoffCol(jend))
403	+	end if
404	+
405	+	groupMaxCutoffCol = 0.0_dp
406	+	gtypeMaxCutoffCol = 0.0_dp
407	+
408	+	#endif
409	+	groupMaxCutoffRow = 0.0_dp
410	+	gtypeMaxCutoffRow = 0.0_dp
411	+
412	+
413	+	!! first we do a single loop over the cutoff groups to find the
414	+	!! largest cutoff for any atypes present in this group. We also
415	+	!! create gtypes at this point.
416	+
417	+	tol = 1.0d-6
418	+	nGroupTypesRow = 0
419	+
420	+	do i = istart, iend
421	+	n_in_i = groupStartRow(i+1) - groupStartRow(i)
422	+	groupMaxCutoffRow(i) = 0.0_dp
423	+	do ia = groupStartRow(i), groupStartRow(i+1)-1
424	+	atom1 = groupListRow(ia)
425	+	#ifdef IS_MPI
426	+	me_i = atid_row(atom1)
427	+	#else
428	+	me_i = atid(atom1)
429	+	#endif
430	+	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
431	+	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
432	+	endif
433	+	enddo
434	+
435	+	if (nGroupTypesRow.eq.0) then
436	+	nGroupTypesRow = nGroupTypesRow + 1
437	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
438	+	groupToGtypeRow(i) = nGroupTypesRow
439	+	else
440	+	GtypeFound = .false.
441	+	do g = 1, nGroupTypesRow
442	+	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
443	+	groupToGtypeRow(i) = g
444	+	GtypeFound = .true.
445	+	endif
446	+	enddo
447	+	if (.not.GtypeFound) then
448	+	nGroupTypesRow = nGroupTypesRow + 1
449	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
450	+	groupToGtypeRow(i) = nGroupTypesRow
451	+	endif
452	+	endif
453	+	enddo
454	+
455	+	#ifdef IS_MPI
456	+	do j = jstart, jend
457	+	n_in_j = groupStartCol(j+1) - groupStartCol(j)
458	+	groupMaxCutoffCol(j) = 0.0_dp
459	+	do ja = groupStartCol(j), groupStartCol(j+1)-1
460	+	atom1 = groupListCol(ja)
461	+
462	+	me_j = atid_col(atom1)
463	+
464	+	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
465	+	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
466	+	endif
467	+	enddo
468	+
469	+	if (nGroupTypesCol.eq.0) then
470	+	nGroupTypesCol = nGroupTypesCol + 1
471	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
472	+	groupToGtypeCol(j) = nGroupTypesCol
473	+	else
474	+	GtypeFound = .false.
475	+	do g = 1, nGroupTypesCol
476	+	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
477	+	groupToGtypeCol(j) = g
478	+	GtypeFound = .true.
479	+	endif
480	+	enddo
481	+	if (.not.GtypeFound) then
482	+	nGroupTypesCol = nGroupTypesCol + 1
483	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
484	+	groupToGtypeCol(j) = nGroupTypesCol
485	+	endif
486	+	endif
487	+	enddo
488	+
489	+	#else
490	+	! Set pointers to information we just found
491	+	nGroupTypesCol = nGroupTypesRow
492	+	groupToGtypeCol => groupToGtypeRow
493	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
494	+	groupMaxCutoffCol => groupMaxCutoffRow
495	+	#endif
496	+
497	+
498	+
499	+
500	+
501	+	!! allocate the gtypeCutoffMap here.
502	+	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
503	+	!! then we do a double loop over all the group TYPES to find the cutoff
504	+	!! map between groups of two types
505	+	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
506	+
507	+	do i = 1, nGroupTypesRow
508	+	do j = 1, nGroupTypesCol
509	+
510	+	select case(cutoffPolicy)
511	+	case(TRADITIONAL_CUTOFF_POLICY)
512	+	thisRcut = tradRcut
513	+	case(MIX_CUTOFF_POLICY)
514	+	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
515	+	case(MAX_CUTOFF_POLICY)
516	+	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
517	+	case default
518	+	call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
519	+	return
520	+	end select
521	+	gtypeCutoffMap(i,j)%rcut = thisRcut
522	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
523	+	skin = defaultRlist - defaultRcut
524	+	listSkin = skin ! set neighbor list skin thickness
525	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
526	+
527	+	! sanity check
528	+
529	+	if (haveDefaultCutoffs) then
530	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
531	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
532	+	endif
533	+	endif
534	+	enddo
535	+	enddo
536	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
537	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
538	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
539	+	#ifdef IS_MPI
540	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
541	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
542	+	#endif
543	+	groupMaxCutoffCol => null()
544	+	gtypeMaxCutoffCol => null()
545	+
546	+	haveGtypeCutoffMap = .true.
547	+	end subroutine createGtypeCutoffMap
548	+
549	+	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
550	+	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
551	+	integer, intent(in) :: cutPolicy
552	+
553	+	defaultRcut = defRcut
554	+	defaultRsw = defRsw
555	+	defaultRlist = defRlist
556	+	cutoffPolicy = cutPolicy
557	+
558	+	haveDefaultCutoffs = .true.
559	+	end subroutine setDefaultCutoffs
560	+
561	+	subroutine setCutoffPolicy(cutPolicy)
562	+
563	+	integer, intent(in) :: cutPolicy
564	+	cutoffPolicy = cutPolicy
565	+	call createGtypeCutoffMap()
566	+	end subroutine setCutoffPolicy
567	+
568	+
569		subroutine setSimVariables()
570		SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
264	–	SIM_uses_LennardJones = SimUsesLennardJones()
265	–	SIM_uses_Electrostatics = SimUsesElectrostatics()
266	–	SIM_uses_Charges = SimUsesCharges()
267	–	SIM_uses_Dipoles = SimUsesDipoles()
268	–	SIM_uses_Sticky = SimUsesSticky()
269	–	SIM_uses_StickyPower = SimUsesStickyPower()
270	–	SIM_uses_GayBerne = SimUsesGayBerne()
571		SIM_uses_EAM = SimUsesEAM()
272	–	SIM_uses_Shapes = SimUsesShapes()
273	–	SIM_uses_FLARB = SimUsesFLARB()
274	–	SIM_uses_RF = SimUsesRF()
572		SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573		SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574		SIM_uses_PBC = SimUsesPBC()
#	Line 288 \| Line 585 \| contains
585
586		error = 0
587
588	<	if (.not. haveInteractionMap) then
588	>	if (.not. haveInteractionHash) then
589	>	myStatus = 0
590	>	call createInteractionHash(myStatus)
591	>	if (myStatus .ne. 0) then
592	>	write(default_error, *) 'createInteractionHash failed in doForces!'
593	>	error = -1
594	>	return
595	>	endif
596	>	endif
597
598	<	myStatus = 0
599	<
600	<	call createInteractionMap(myStatus)
296	<
598	>	if (.not. haveGtypeCutoffMap) then
599	>	myStatus = 0
600	>	call createGtypeCutoffMap(myStatus)
601		if (myStatus .ne. 0) then
602	<	write(default_error, *) 'createInteractionMap failed in doForces!'
602	>	write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
603		error = -1
604		return
605		endif
#	Line 305 \| Line 609 \| contains
609		call setSimVariables()
610		endif
611
612	<	if (.not. haveRlist) then
613	<	write(default_error, *) 'rList has not been set in doForces!'
614	<	error = -1
615	<	return
616	<	endif
612	>	! if (.not. haveRlist) then
613	>	! write(default_error, *) 'rList has not been set in doForces!'
614	>	! error = -1
615	>	! return
616	>	! endif
617
618		if (.not. haveNeighborList) then
619		write(default_error, *) 'neighbor list has not been initialized in doForces!'
#	Line 334 \| Line 638 \| contains
638		end subroutine doReadyCheck
639
640
641	<	subroutine init_FF(use_RF_c, thisStat)
641	>	subroutine init_FF(thisESM, thisStat)
642
643	<	logical, intent(in) :: use_RF_c
340	<
643	>	integer, intent(in) :: thisESM
644		integer, intent(out) :: thisStat
645		integer :: my_status, nMatches
646		integer, pointer :: MatchList(:) => null()
344	–	real(kind=dp) :: rcut, rrf, rt, dielect
647
648		!! assume things are copacetic, unless they aren't
649		thisStat = 0
650
651	<	!! Fortran's version of a cast:
350	<	FF_uses_RF = use_RF_c
651	>	electrostaticSummationMethod = thisESM
652
653		!! init_FF is called after all of the atom types have been
654		!! defined in atype_module using the new_atype subroutine.
#	Line 356 \| Line 657 \| contains
657		!! interactions are used by the force field.
658
659		FF_uses_DirectionalAtoms = .false.
359	–	FF_uses_LennardJones = .false.
360	–	FF_uses_Electrostatics = .false.
361	–	FF_uses_Charges = .false.
660		FF_uses_Dipoles = .false.
363	–	FF_uses_Sticky = .false.
364	–	FF_uses_StickyPower = .false.
661		FF_uses_GayBerne = .false.
662		FF_uses_EAM = .false.
367	–	FF_uses_Shapes = .false.
368	–	FF_uses_FLARB = .false.
663
664		call getMatchingElementList(atypes, "is_Directional", .true., &
665		nMatches, MatchList)
666		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
667
374	–	call getMatchingElementList(atypes, "is_LennardJones", .true., &
375	–	nMatches, MatchList)
376	–	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
377	–
378	–	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
379	–	nMatches, MatchList)
380	–	if (nMatches .gt. 0) then
381	–	FF_uses_Electrostatics = .true.
382	–	endif
383	–
384	–	call getMatchingElementList(atypes, "is_Charge", .true., &
385	–	nMatches, MatchList)
386	–	if (nMatches .gt. 0) then
387	–	FF_uses_Charges = .true.
388	–	FF_uses_Electrostatics = .true.
389	–	endif
390	–
668		call getMatchingElementList(atypes, "is_Dipole", .true., &
669		nMatches, MatchList)
670	<	if (nMatches .gt. 0) then
394	<	FF_uses_Dipoles = .true.
395	<	FF_uses_Electrostatics = .true.
396	<	FF_uses_DirectionalAtoms = .true.
397	<	endif
398	<
399	<	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
400	<	nMatches, MatchList)
401	<	if (nMatches .gt. 0) then
402	<	FF_uses_Quadrupoles = .true.
403	<	FF_uses_Electrostatics = .true.
404	<	FF_uses_DirectionalAtoms = .true.
405	<	endif
406	<
407	<	call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
408	<	MatchList)
409	<	if (nMatches .gt. 0) then
410	<	FF_uses_Sticky = .true.
411	<	FF_uses_DirectionalAtoms = .true.
412	<	endif
413	<
414	<	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
415	<	MatchList)
416	<	if (nMatches .gt. 0) then
417	<	FF_uses_StickyPower = .true.
418	<	FF_uses_DirectionalAtoms = .true.
419	<	endif
670	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
671
672		call getMatchingElementList(atypes, "is_GayBerne", .true., &
673		nMatches, MatchList)
674	<	if (nMatches .gt. 0) then
424	<	FF_uses_GayBerne = .true.
425	<	FF_uses_DirectionalAtoms = .true.
426	<	endif
674	>	if (nMatches .gt. 0) FF_uses_GayBerne = .true.
675
676		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
677		if (nMatches .gt. 0) FF_uses_EAM = .true.
430	–
431	–	call getMatchingElementList(atypes, "is_Shape", .true., &
432	–	nMatches, MatchList)
433	–	if (nMatches .gt. 0) then
434	–	FF_uses_Shapes = .true.
435	–	FF_uses_DirectionalAtoms = .true.
436	–	endif
678
438	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
439	–	nMatches, MatchList)
440	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
679
442	–	!! Assume sanity (for the sake of argument)
680		haveSaneForceField = .true.
681
445	–	!! check to make sure the FF_uses_RF setting makes sense
446	–
447	–	if (FF_uses_dipoles) then
448	–	if (FF_uses_RF) then
449	–	dielect = getDielect()
450	–	call initialize_rf(dielect)
451	–	endif
452	–	else
453	–	if (FF_uses_RF) then
454	–	write(default_error,*) 'Using Reaction Field with no dipoles? Huh?'
455	–	thisStat = -1
456	–	haveSaneForceField = .false.
457	–	return
458	–	endif
459	–	endif
460	–
461	–	!sticky module does not contain check_sticky_FF anymore
462	–	!if (FF_uses_sticky) then
463	–	! call check_sticky_FF(my_status)
464	–	! if (my_status /= 0) then
465	–	! thisStat = -1
466	–	! haveSaneForceField = .false.
467	–	! return
468	–	! end if
469	–	!endif
470	–
682		if (FF_uses_EAM) then
683		call init_EAM_FF(my_status)
684		if (my_status /= 0) then
#	Line 478 \| Line 689 \| contains
689		end if
690		endif
691
481	–	if (FF_uses_GayBerne) then
482	–	call check_gb_pair_FF(my_status)
483	–	if (my_status .ne. 0) then
484	–	thisStat = -1
485	–	haveSaneForceField = .false.
486	–	return
487	–	endif
488	–	endif
489	–
490	–	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
491	–	endif
492	–
692		if (.not. haveNeighborList) then
693		!! Create neighbor lists
694		call expandNeighborList(nLocal, my_status)
#	Line 523 \| Line 722 \| contains
722
723		!! Stress Tensor
724		real( kind = dp), dimension(9) :: tau
725	<	real ( kind = dp ) :: pot
725	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
726		logical ( kind = 2) :: do_pot_c, do_stress_c
727		logical :: do_pot
728		logical :: do_stress
729		logical :: in_switching_region
730		#ifdef IS_MPI
731	<	real( kind = DP ) :: pot_local
731	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
732		integer :: nAtomsInRow
733		integer :: nAtomsInCol
734		integer :: nprocs
#	Line 544 \| Line 743 \| contains
743		integer :: nlist
744		real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
745		real( kind = DP ) :: sw, dswdr, swderiv, mf
746	+	real( kind = DP ) :: rVal
747		real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
748		real(kind=dp) :: rfpot, mu_i, virial
749		integer :: me_i, me_j, n_in_i, n_in_j
#	Line 553 \| Line 753 \| contains
753		integer :: localError
754		integer :: propPack_i, propPack_j
755		integer :: loopStart, loopEnd, loop
756	+	integer :: iHash
757	+	integer :: i1
758	+
759
557	–	real(kind=dp) :: listSkin = 1.0
558	–
760		!! initialize local variables
761
762		#ifdef IS_MPI
#	Line 672 \| Line 873 \| contains
873		endif
874
875		#ifdef IS_MPI
876	+	me_j = atid_col(j)
877		call get_interatomic_vector(q_group_Row(:,i), &
878		q_group_Col(:,j), d_grp, rgrpsq)
879		#else
880	+	me_j = atid(j)
881		call get_interatomic_vector(q_group(:,i), &
882		q_group(:,j), d_grp, rgrpsq)
883	<	#endif
883	>	#endif
884
885	<	if (rgrpsq < rlistsq) then
885	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
886		if (update_nlist) then
887		nlist = nlist + 1
888
#	Line 699 \| Line 902 \| contains
902
903		list(nlist) = j
904		endif
905	<
906	<	if (loop .eq. PAIR_LOOP) then
704	<	vij = 0.0d0
705	<	fij(1:3) = 0.0d0
706	<	endif
707	<
708	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
709	<	in_switching_region)
710	<
711	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
712	<
713	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
714	<
715	<	atom1 = groupListRow(ia)
716	<
717	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
905	>
906	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
907
908	<	atom2 = groupListCol(jb)
909	<
910	<	if (skipThisPair(atom1, atom2)) cycle inner
911	<
912	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
913	<	d_atm(1:3) = d_grp(1:3)
914	<	ratmsq = rgrpsq
915	<	else
908	>	if (loop .eq. PAIR_LOOP) then
909	>	vij = 0.0d0
910	>	fij(1:3) = 0.0d0
911	>	endif
912	>
913	>	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
914	>	in_switching_region)
915	>
916	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
917	>
918	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
919	>
920	>	atom1 = groupListRow(ia)
921	>
922	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
923	>
924	>	atom2 = groupListCol(jb)
925	>
926	>	if (skipThisPair(atom1, atom2)) cycle inner
927	>
928	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
929	>	d_atm(1:3) = d_grp(1:3)
930	>	ratmsq = rgrpsq
931	>	else
932		#ifdef IS_MPI
933	<	call get_interatomic_vector(q_Row(:,atom1), &
934	<	q_Col(:,atom2), d_atm, ratmsq)
933	>	call get_interatomic_vector(q_Row(:,atom1), &
934	>	q_Col(:,atom2), d_atm, ratmsq)
935		#else
936	<	call get_interatomic_vector(q(:,atom1), &
937	<	q(:,atom2), d_atm, ratmsq)
936	>	call get_interatomic_vector(q(:,atom1), &
937	>	q(:,atom2), d_atm, ratmsq)
938		#endif
939	<	endif
940	<
941	<	if (loop .eq. PREPAIR_LOOP) then
939	>	endif
940	>
941	>	if (loop .eq. PREPAIR_LOOP) then
942		#ifdef IS_MPI
943	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
944	<	rgrpsq, d_grp, do_pot, do_stress, &
945	<	eFrame, A, f, t, pot_local)
943	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
944	>	rgrpsq, d_grp, do_pot, do_stress, &
945	>	eFrame, A, f, t, pot_local)
946		#else
947	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
948	<	rgrpsq, d_grp, do_pot, do_stress, &
949	<	eFrame, A, f, t, pot)
947	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
948	>	rgrpsq, d_grp, do_pot, do_stress, &
949	>	eFrame, A, f, t, pot)
950		#endif
951	<	else
951	>	else
952		#ifdef IS_MPI
953	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
954	<	do_pot, &
955	<	eFrame, A, f, t, pot_local, vpair, fpair)
953	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
954	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
955	>	fpair, d_grp, rgrp)
956		#else
957	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
958	<	do_pot, &
959	<	eFrame, A, f, t, pot, vpair, fpair)
957	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
958	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
959	>	d_grp, rgrp)
960		#endif
961	+	vij = vij + vpair
962	+	fij(1:3) = fij(1:3) + fpair(1:3)
963	+	endif
964	+	enddo inner
965	+	enddo
966
967	<	vij = vij + vpair
968	<	fij(1:3) = fij(1:3) + fpair(1:3)
969	<	endif
970	<	enddo inner
971	<	enddo
972	<
973	<	if (loop .eq. PAIR_LOOP) then
974	<	if (in_switching_region) then
975	<	swderiv = vij*dswdr/rgrp
976	<	fij(1) = fij(1) + swderiv*d_grp(1)
767	<	fij(2) = fij(2) + swderiv*d_grp(2)
768	<	fij(3) = fij(3) + swderiv*d_grp(3)
769	<
770	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
771	<	atom1=groupListRow(ia)
772	<	mf = mfactRow(atom1)
967	>	if (loop .eq. PAIR_LOOP) then
968	>	if (in_switching_region) then
969	>	swderiv = vij*dswdr/rgrp
970	>	fij(1) = fij(1) + swderiv*d_grp(1)
971	>	fij(2) = fij(2) + swderiv*d_grp(2)
972	>	fij(3) = fij(3) + swderiv*d_grp(3)
973	>
974	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
975	>	atom1=groupListRow(ia)
976	>	mf = mfactRow(atom1)
977		#ifdef IS_MPI
978	<	f_Row(1,atom1) = f_Row(1,atom1) + swderivd_grp(1)mf
979	<	f_Row(2,atom1) = f_Row(2,atom1) + swderivd_grp(2)mf
980	<	f_Row(3,atom1) = f_Row(3,atom1) + swderivd_grp(3)mf
978	>	f_Row(1,atom1) = f_Row(1,atom1) + swderivd_grp(1)mf
979	>	f_Row(2,atom1) = f_Row(2,atom1) + swderivd_grp(2)mf
980	>	f_Row(3,atom1) = f_Row(3,atom1) + swderivd_grp(3)mf
981		#else
982	<	f(1,atom1) = f(1,atom1) + swderivd_grp(1)mf
983	<	f(2,atom1) = f(2,atom1) + swderivd_grp(2)mf
984	<	f(3,atom1) = f(3,atom1) + swderivd_grp(3)mf
982	>	f(1,atom1) = f(1,atom1) + swderivd_grp(1)mf
983	>	f(2,atom1) = f(2,atom1) + swderivd_grp(2)mf
984	>	f(3,atom1) = f(3,atom1) + swderivd_grp(3)mf
985		#endif
986	<	enddo
987	<
988	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
989	<	atom2=groupListCol(jb)
990	<	mf = mfactCol(atom2)
986	>	enddo
987	>
988	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
989	>	atom2=groupListCol(jb)
990	>	mf = mfactCol(atom2)
991		#ifdef IS_MPI
992	<	f_Col(1,atom2) = f_Col(1,atom2) - swderivd_grp(1)mf
993	<	f_Col(2,atom2) = f_Col(2,atom2) - swderivd_grp(2)mf
994	<	f_Col(3,atom2) = f_Col(3,atom2) - swderivd_grp(3)mf
992	>	f_Col(1,atom2) = f_Col(1,atom2) - swderivd_grp(1)mf
993	>	f_Col(2,atom2) = f_Col(2,atom2) - swderivd_grp(2)mf
994	>	f_Col(3,atom2) = f_Col(3,atom2) - swderivd_grp(3)mf
995		#else
996	<	f(1,atom2) = f(1,atom2) - swderivd_grp(1)mf
997	<	f(2,atom2) = f(2,atom2) - swderivd_grp(2)mf
998	<	f(3,atom2) = f(3,atom2) - swderivd_grp(3)mf
996	>	f(1,atom2) = f(1,atom2) - swderivd_grp(1)mf
997	>	f(2,atom2) = f(2,atom2) - swderivd_grp(2)mf
998	>	f(3,atom2) = f(3,atom2) - swderivd_grp(3)mf
999		#endif
1000	<	enddo
1001	<	endif
1000	>	enddo
1001	>	endif
1002
1003	<	if (do_stress) call add_stress_tensor(d_grp, fij)
1003	>	if (do_stress) call add_stress_tensor(d_grp, fij)
1004	>	endif
1005		endif
1006	<	end if
1006	>	endif
1007		enddo
1008	+
1009		enddo outer
1010
1011		if (update_nlist) then
#	Line 859 \| Line 1065 \| contains
1065
1066		if (do_pot) then
1067		! scatter/gather pot_row into the members of my column
1068	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1069	<
1068	>	do i = 1,LR_POT_TYPES
1069	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1070	>	end do
1071		! scatter/gather pot_local into all other procs
1072		! add resultant to get total pot
1073		do i = 1, nlocal
1074	<	pot_local = pot_local + pot_Temp(i)
1074	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1075	>	+ pot_Temp(1:LR_POT_TYPES,i)
1076		enddo
1077
1078		pot_Temp = 0.0_DP
1079	<
1080	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1079	>	do i = 1,LR_POT_TYPES
1080	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1081	>	end do
1082		do i = 1, nlocal
1083	<	pot_local = pot_local + pot_Temp(i)
1083	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1084	>	+ pot_Temp(1:LR_POT_TYPES,i)
1085		enddo
1086
1087		endif
1088		#endif
1089
1090	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1091	<
1092	<	if (FF_uses_RF .and. SIM_uses_RF) then
1090	>	if (SIM_requires_postpair_calc) then
1091	>	do i = 1, nlocal
1092	>
1093	>	! we loop only over the local atoms, so we don't need row and column
1094	>	! lookups for the types
1095	>
1096	>	me_i = atid(i)
1097	>
1098	>	! is the atom electrostatic? See if it would have an
1099	>	! electrostatic interaction with itself
1100	>	iHash = InteractionHash(me_i,me_i)
1101
1102	+	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1103		#ifdef IS_MPI
1104	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1105	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
887	<	do i = 1,nlocal
888	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
889	<	end do
890	<	#endif
891	<
892	<	do i = 1, nLocal
893	<
894	<	rfpot = 0.0_DP
895	<	#ifdef IS_MPI
896	<	me_i = atid_row(i)
1104	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1105	>	t, do_pot)
1106		#else
1107	<	me_i = atid(i)
1107	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1108	>	t, do_pot)
1109		#endif
1110	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1110	>	endif
1111	>
1112	>
1113	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1114
1115	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1116	<
1117	<	mu_i = getDipoleMoment(me_i)
1118	<
1119	<	!! The reaction field needs to include a self contribution
1120	<	!! to the field:
1121	<	call accumulate_self_rf(i, mu_i, eFrame)
1122	<	!! Get the reaction field contribution to the
1123	<	!! potential and torques:
1124	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1115	>	! loop over the excludes to accumulate RF stuff we've
1116	>	! left out of the normal pair loop
1117	>
1118	>	do i1 = 1, nSkipsForAtom(i)
1119	>	j = skipsForAtom(i, i1)
1120	>
1121	>	! prevent overcounting of the skips
1122	>	if (i.lt.j) then
1123	>	call get_interatomic_vector(q(:,i), &
1124	>	q(:,j), d_atm, ratmsq)
1125	>	rVal = dsqrt(ratmsq)
1126	>	call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1127	>	in_switching_region)
1128		#ifdef IS_MPI
1129	<	pot_local = pot_local + rfpot
1129	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1130	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1131		#else
1132	<	pot = pot + rfpot
1133	<
1132	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1133	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1134		#endif
1135	<	endif
1136	<	enddo
1137	<	endif
1135	>	endif
1136	>	enddo
1137	>	endif
1138	>	enddo
1139		endif
1140	<
923	<
1140	>
1141		#ifdef IS_MPI
1142	<
1142	>
1143		if (do_pot) then
1144	<	pot = pot + pot_local
1145	<	!! we assume the c code will do the allreduce to get the total potential
929	<	!! we could do it right here if we needed to...
1144	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1145	>	mpi_comm_world,mpi_err)
1146		endif
1147	<
1147	>
1148		if (do_stress) then
1149		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1150		mpi_comm_world,mpi_err)
1151		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1152		mpi_comm_world,mpi_err)
1153		endif
1154	<
1154	>
1155		#else
1156	<
1156	>
1157		if (do_stress) then
1158		tau = tau_Temp
1159		virial = virial_Temp
1160		endif
1161	<
1161	>
1162		#endif
1163	<
1163	>
1164		end subroutine do_force_loop
1165
1166		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1167	<	eFrame, A, f, t, pot, vpair, fpair)
1167	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1168
1169	<	real( kind = dp ) :: pot, vpair, sw
1169	>	real( kind = dp ) :: vpair, sw
1170	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1171		real( kind = dp ), dimension(3) :: fpair
1172		real( kind = dp ), dimension(nLocal) :: mfact
1173		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 961 \| Line 1178 \| contains
1178		logical, intent(inout) :: do_pot
1179		integer, intent(in) :: i, j
1180		real ( kind = dp ), intent(inout) :: rijsq
1181	<	real ( kind = dp ) :: r
1181	>	real ( kind = dp ), intent(inout) :: r_grp
1182		real ( kind = dp ), intent(inout) :: d(3)
1183	<	real ( kind = dp ) :: ebalance
1183	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1184	>	real ( kind = dp ) :: r
1185		integer :: me_i, me_j
1186
1187	<	integer :: iMap
1187	>	integer :: iHash
1188
1189		r = sqrt(rijsq)
1190		vpair = 0.0d0
#	Line 980 \| Line 1198 \| contains
1198		me_j = atid(j)
1199		#endif
1200
1201	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1202	<
1203	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1204	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1201	>	iHash = InteractionHash(me_i, me_j)
1202	>
1203	>	if ( iand(iHash, LJ_PAIR).ne.0 ) then
1204	>	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1205	>	pot(VDW_POT), f, do_pot)
1206		endif
1207	<
1208	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1207	>
1208	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1209		call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1210	<	pot, eFrame, f, t, do_pot)
992	<
993	<	if (FF_uses_RF .and. SIM_uses_RF) then
994	<
995	<	! CHECK ME (RF needs to know about all electrostatic types)
996	<	call accumulate_rf(i, j, r, eFrame, sw)
997	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
998	<	endif
999	<
1210	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1211		endif
1212	<
1213	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1212	>
1213	>	if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1214		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1215	<	pot, A, f, t, do_pot)
1215	>	pot(HB_POT), A, f, t, do_pot)
1216		endif
1217	<
1218	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1217	>
1218	>	if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1219		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1220	<	pot, A, f, t, do_pot)
1220	>	pot(HB_POT), A, f, t, do_pot)
1221		endif
1222	<
1223	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1222	>
1223	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1224		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1225	<	pot, A, f, t, do_pot)
1225	>	pot(VDW_POT), A, f, t, do_pot)
1226		endif
1227
1228	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1229	<	call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1230	<	pot, A, f, t, do_pot)
1228	>	if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1229	>	call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1230	>	pot(VDW_POT), A, f, t, do_pot)
1231		endif
1232	<
1233	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1234	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1235	<	do_pot)
1232	>
1233	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1234	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1235	>	pot(METALLIC_POT), f, do_pot)
1236		endif
1237	<
1238	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1237	>
1238	>	if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1239		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1240	<	pot, A, f, t, do_pot)
1240	>	pot(VDW_POT), A, f, t, do_pot)
1241		endif
1242	<
1243	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1242	>
1243	>	if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1244		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1245	<	pot, A, f, t, do_pot)
1245	>	pot(VDW_POT), A, f, t, do_pot)
1246		endif
1247	<
1247	>
1248		end subroutine do_pair
1249
1250		subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1251		do_pot, do_stress, eFrame, A, f, t, pot)
1252
1253	<	real( kind = dp ) :: pot, sw
1253	>	real( kind = dp ) :: sw
1254	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1255		real( kind = dp ), dimension(9,nLocal) :: eFrame
1256		real (kind=dp), dimension(9,nLocal) :: A
1257		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1051 \| Line 1263 \| contains
1263		real ( kind = dp ) :: r, rc
1264		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1265
1266	<	integer :: me_i, me_j, iMap
1266	>	integer :: me_i, me_j, iHash
1267
1268	+	r = sqrt(rijsq)
1269	+
1270		#ifdef IS_MPI
1271		me_i = atid_row(i)
1272		me_j = atid_col(j)
#	Line 1061 \| Line 1275 \| contains
1275		me_j = atid(j)
1276		#endif
1277
1278	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1278	>	iHash = InteractionHash(me_i, me_j)
1279
1280	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1280	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1281		call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1282		endif
1283
#	Line 1072 \| Line 1286 \| contains
1286
1287		subroutine do_preforce(nlocal,pot)
1288		integer :: nlocal
1289	<	real( kind = dp ) :: pot
1289	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1290
1291		if (FF_uses_EAM .and. SIM_uses_EAM) then
1292	<	call calc_EAM_preforce_Frho(nlocal,pot)
1292	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1293		endif
1294
1295
#	Line 1161 \| Line 1375 \| contains
1375		pot_Col = 0.0_dp
1376		pot_Temp = 0.0_dp
1377
1164	–	rf_Row = 0.0_dp
1165	–	rf_Col = 0.0_dp
1166	–	rf_Temp = 0.0_dp
1167	–
1378		#endif
1379
1380		if (FF_uses_EAM .and. SIM_uses_EAM) then
1381		call clean_EAM()
1382		endif
1383
1174	–	rf = 0.0_dp
1384		tau_Temp = 0.0_dp
1385		virial_Temp = 0.0_dp
1386		end subroutine zero_work_arrays
#	Line 1260 \| Line 1469 \| contains
1469
1470		function FF_UsesDirectionalAtoms() result(doesit)
1471		logical :: doesit
1472	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1264	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1265	<	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1472	>	doesit = FF_uses_DirectionalAtoms
1473		end function FF_UsesDirectionalAtoms
1474
1475		function FF_RequiresPrepairCalc() result(doesit)
#	Line 1270 \| Line 1477 \| contains
1477		doesit = FF_uses_EAM
1478		end function FF_RequiresPrepairCalc
1479
1273	–	function FF_RequiresPostpairCalc() result(doesit)
1274	–	logical :: doesit
1275	–	doesit = FF_uses_RF
1276	–	end function FF_RequiresPostpairCalc
1277	–
1480		#ifdef PROFILE
1481		function getforcetime() result(totalforcetime)
1482		real(kind=dp) :: totalforcetime

Diff Legend

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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents): Revision 2261 by gezelter, Tue Jun 28 13:58:45 2005 UTC vs. Revision 2411 by chrisfen, Wed Nov 2 21:01:21 2005 UTC

Diff Legend

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2261 by gezelter, Tue Jun 28 13:58:45 2005 UTC vs.
Revision 2411 by chrisfen, Wed Nov 2 21:01:21 2005 UTC