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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2270 by gezelter, Tue Aug 9 22:33:37 2005 UTC vs.
Revision 2342 by chrisfen, Tue Oct 4 19:32:58 2005 UTC

#	Line 45 \| Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.28 2005-08-09 22:33:37 gezelter Exp $, $Date: 2005-08-09 22:33:37 $, $Name: not supported by cvs2svn $, $Revision: 1.28 $
48	>	!! @version $Id: doForces.F90,v 1.50 2005-10-04 19:32:58 chrisfen Exp $, $Date: 2005-10-04 19:32:58 $, $Name: not supported by cvs2svn $, $Revision: 1.50 $
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
61	>	use reaction_field_module
62		use gb_pair
63		use shapes
64		use vector_class
#	Line 73 \| Line 73 \| module doForces
73
74		#define __FORTRAN90
75		#include "UseTheForce/fSwitchingFunction.h"
76	+	#include "UseTheForce/fCutoffPolicy.h"
77		#include "UseTheForce/DarkSide/fInteractionMap.h"
78	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79
80	+
81		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
82		INTEGER, PARAMETER:: PAIR_LOOP = 2
83
#	Line 83 \| Line 86 \| module doForces
86		logical, save :: haveSaneForceField = .false.
87		logical, save :: haveInteractionHash = .false.
88		logical, save :: haveGtypeCutoffMap = .false.
89	+	logical, save :: haveDefaultCutoffs = .false.
90	+	logical, save :: haveRlist = .false.
91
92		logical, save :: FF_uses_DirectionalAtoms
93		logical, save :: FF_uses_Dipoles
94		logical, save :: FF_uses_GayBerne
95		logical, save :: FF_uses_EAM
91	–	logical, save :: FF_uses_RF
96
97		logical, save :: SIM_uses_DirectionalAtoms
98		logical, save :: SIM_uses_EAM
95	–	logical, save :: SIM_uses_RF
99		logical, save :: SIM_requires_postpair_calc
100		logical, save :: SIM_requires_prepair_calc
101		logical, save :: SIM_uses_PBC
102
103	+	integer, save :: electrostaticSummationMethod
104	+
105		public :: init_FF
106	+	public :: setDefaultCutoffs
107		public :: do_force_loop
108		public :: createInteractionHash
109		public :: createGtypeCutoffMap
110	+	public :: getStickyCut
111	+	public :: getStickyPowerCut
112	+	public :: getGayBerneCut
113	+	public :: getEAMCut
114	+	public :: getShapeCut
115
116		#ifdef PROFILE
117		public :: getforcetime
#	Line 122 \| Line 133 \| module doForces
133		real(kind=dp) :: rlistsq
134		end type gtypeCutoffs
135		type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
136	+
137	+	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
138	+	real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
139	+	real(kind=dp),save :: listSkin
140
141		contains
142
#	Line 146 \| Line 161 \| contains
161		logical :: j_is_GB
162		logical :: j_is_EAM
163		logical :: j_is_Shape
164	<
164	>	real(kind=dp) :: myRcut
165	>
166		status = 0
167
168		if (.not. associated(atypes)) then
#	Line 179 \| Line 195 \| contains
195		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
196		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
197
182	–	if (i_is_LJ) then
183	–	thisCut = getDefaultLJCutoff(i)
184	–	if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
185	–	endif
186	–
187	–
188	–
198		do j = i, nAtypes
199
200		iHash = 0
#	Line 238 \| Line 247 \| contains
247		haveInteractionHash = .true.
248		end subroutine createInteractionHash
249
250	<	subroutine createGtypeCutoffMap(defaultRcut, defaultSkinThickness, stat)
242	<
243	<	real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
244	<	integer, intent(out) :: stat
250	>	subroutine createGtypeCutoffMap(stat)
251
252	<	integer :: myStatus, nAtypes
252	>	integer, intent(out), optional :: stat
253	>	logical :: i_is_LJ
254	>	logical :: i_is_Elect
255	>	logical :: i_is_Sticky
256	>	logical :: i_is_StickyP
257	>	logical :: i_is_GB
258	>	logical :: i_is_EAM
259	>	logical :: i_is_Shape
260	>	logical :: GtypeFound
261
262	+	integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend
263	+	integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
264	+	integer :: nGroupsInRow
265	+	real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
266	+	real(kind=dp) :: biggestAtypeCutoff
267	+
268		stat = 0
269		if (.not. haveInteractionHash) then
270		call createInteractionHash(myStatus)
#	Line 254 \| Line 274 \| contains
274		return
275		endif
276		endif
277	<
277	>	#ifdef IS_MPI
278	>	nGroupsInRow = getNgroupsInRow(plan_group_row)
279	>	#endif
280		nAtypes = getSize(atypes)
281	<
281	>	! Set all of the initial cutoffs to zero.
282	>	atypeMaxCutoff = 0.0_dp
283		do i = 1, nAtypes
284	<
285	<	atypeMaxCutoff(i) =
284	>	if (SimHasAtype(i)) then
285	>	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
286	>	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
287	>	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
288	>	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
289	>	call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
290	>	call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
291	>	call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
292	>
293	>
294	>	if (haveDefaultCutoffs) then
295	>	atypeMaxCutoff(i) = defaultRcut
296	>	else
297	>	if (i_is_LJ) then
298	>	thisRcut = getSigma(i) * 2.5_dp
299	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
300	>	endif
301	>	if (i_is_Elect) then
302	>	thisRcut = defaultRcut
303	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
304	>	endif
305	>	if (i_is_Sticky) then
306	>	thisRcut = getStickyCut(i)
307	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308	>	endif
309	>	if (i_is_StickyP) then
310	>	thisRcut = getStickyPowerCut(i)
311	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312	>	endif
313	>	if (i_is_GB) then
314	>	thisRcut = getGayBerneCut(i)
315	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316	>	endif
317	>	if (i_is_EAM) then
318	>	thisRcut = getEAMCut(i)
319	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
320	>	endif
321	>	if (i_is_Shape) then
322	>	thisRcut = getShapeCut(i)
323	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
324	>	endif
325	>	endif
326	>
327	>
328	>	if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
329	>	biggestAtypeCutoff = atypeMaxCutoff(i)
330	>	endif
331
332	+	endif
333	+	enddo
334	+
335	+	nGroupTypes = 0
336
337	+	istart = 1
338	+	#ifdef IS_MPI
339	+	iend = nGroupsInRow
340	+	#else
341	+	iend = nGroups
342	+	#endif
343	+
344	+	!! allocate the groupToGtype and gtypeMaxCutoff here.
345	+	if(.not.allocated(groupToGtype)) then
346	+	allocate(groupToGtype(iend))
347	+	allocate(groupMaxCutoff(iend))
348	+	allocate(gtypeMaxCutoff(iend))
349	+	groupMaxCutoff = 0.0_dp
350	+	gtypeMaxCutoff = 0.0_dp
351	+	endif
352	+	!! first we do a single loop over the cutoff groups to find the
353	+	!! largest cutoff for any atypes present in this group. We also
354	+	!! create gtypes at this point.
355	+
356	+	tol = 1.0d-6
357	+
358	+	do i = istart, iend
359	+	n_in_i = groupStartRow(i+1) - groupStartRow(i)
360	+	groupMaxCutoff(i) = 0.0_dp
361	+	do ia = groupStartRow(i), groupStartRow(i+1)-1
362	+	atom1 = groupListRow(ia)
363	+	#ifdef IS_MPI
364	+	me_i = atid_row(atom1)
365	+	#else
366	+	me_i = atid(atom1)
367	+	#endif
368	+	if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
369	+	groupMaxCutoff(i)=atypeMaxCutoff(me_i)
370	+	endif
371	+	enddo
372
373	+	if (nGroupTypes.eq.0) then
374	+	nGroupTypes = nGroupTypes + 1
375	+	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
376	+	groupToGtype(i) = nGroupTypes
377	+	else
378	+	GtypeFound = .false.
379	+	do g = 1, nGroupTypes
380	+	if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
381	+	groupToGtype(i) = g
382	+	GtypeFound = .true.
383	+	endif
384	+	enddo
385	+	if (.not.GtypeFound) then
386	+	nGroupTypes = nGroupTypes + 1
387	+	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
388	+	groupToGtype(i) = nGroupTypes
389	+	endif
390	+	endif
391	+	enddo
392
393	+	!! allocate the gtypeCutoffMap here.
394	+	allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
395	+	!! then we do a double loop over all the group TYPES to find the cutoff
396	+	!! map between groups of two types
397	+
398	+	do i = 1, nGroupTypes
399	+	do j = 1, nGroupTypes
400	+
401	+	select case(cutoffPolicy)
402	+	case(TRADITIONAL_CUTOFF_POLICY)
403	+	thisRcut = maxval(gtypeMaxCutoff)
404	+	case(MIX_CUTOFF_POLICY)
405	+	thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
406	+	case(MAX_CUTOFF_POLICY)
407	+	thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
408	+	case default
409	+	call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
410	+	return
411	+	end select
412	+	gtypeCutoffMap(i,j)%rcut = thisRcut
413	+	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
414	+	skin = defaultRlist - defaultRcut
415	+	listSkin = skin ! set neighbor list skin thickness
416	+	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
417
418	<	haveGtypeCutoffMap = .true.
418	>	! sanity check
419	>
420	>	if (haveDefaultCutoffs) then
421	>	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
422	>	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
423	>	endif
424	>	endif
425	>	enddo
426	>	enddo
427	>
428	>	haveGtypeCutoffMap = .true.
429		end subroutine createGtypeCutoffMap
430	+
431	+	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
432	+	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
433	+	integer, intent(in) :: cutPolicy
434
435	+	defaultRcut = defRcut
436	+	defaultRsw = defRsw
437	+	defaultRlist = defRlist
438	+	cutoffPolicy = cutPolicy
439	+
440	+	haveDefaultCutoffs = .true.
441	+	end subroutine setDefaultCutoffs
442	+
443	+	subroutine setCutoffPolicy(cutPolicy)
444	+
445	+	integer, intent(in) :: cutPolicy
446	+	cutoffPolicy = cutPolicy
447	+	call createGtypeCutoffMap()
448	+	end subroutine setCutoffPolicy
449	+
450	+
451		subroutine setSimVariables()
452		SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
453		SIM_uses_EAM = SimUsesEAM()
274	–	SIM_uses_RF = SimUsesRF()
454		SIM_requires_postpair_calc = SimRequiresPostpairCalc()
455		SIM_requires_prepair_calc = SimRequiresPrepairCalc()
456		SIM_uses_PBC = SimUsesPBC()
#	Line 312 \| Line 491 \| contains
491		call setSimVariables()
492		endif
493
494	<	if (.not. haveRlist) then
495	<	write(default_error, *) 'rList has not been set in doForces!'
496	<	error = -1
497	<	return
498	<	endif
494	>	! if (.not. haveRlist) then
495	>	! write(default_error, *) 'rList has not been set in doForces!'
496	>	! error = -1
497	>	! return
498	>	! endif
499
500		if (.not. haveNeighborList) then
501		write(default_error, *) 'neighbor list has not been initialized in doForces!'
#	Line 341 \| Line 520 \| contains
520		end subroutine doReadyCheck
521
522
523	<	subroutine init_FF(use_RF_c, thisStat)
523	>	subroutine init_FF(thisESM, thisStat)
524
525	<	logical, intent(in) :: use_RF_c
347	<
525	>	integer, intent(in) :: thisESM
526		integer, intent(out) :: thisStat
527		integer :: my_status, nMatches
528		integer, pointer :: MatchList(:) => null()
#	Line 353 \| Line 531 \| contains
531		!! assume things are copacetic, unless they aren't
532		thisStat = 0
533
534	<	!! Fortran's version of a cast:
357	<	FF_uses_RF = use_RF_c
534	>	electrostaticSummationMethod = thisESM
535
536		!! init_FF is called after all of the atom types have been
537		!! defined in atype_module using the new_atype subroutine.
#	Line 385 \| Line 562 \| contains
562
563		haveSaneForceField = .true.
564
565	<	!! check to make sure the FF_uses_RF setting makes sense
565	>	!! check to make sure the reaction field setting makes sense
566
567		if (FF_uses_Dipoles) then
568	<	if (FF_uses_RF) then
568	>	if (electrostaticSummationMethod == REACTION_FIELD) then
569		dielect = getDielect()
570		call initialize_rf(dielect)
571		endif
572		else
573	<	if (FF_uses_RF) then
573	>	if (electrostaticSummationMethod == REACTION_FIELD) then
574		write(default_error,*) 'Using Reaction Field with no dipoles? Huh?'
575		thisStat = -1
576		haveSaneForceField = .false.
#	Line 484 \| Line 661 \| contains
661		integer :: propPack_i, propPack_j
662		integer :: loopStart, loopEnd, loop
663		integer :: iHash
664	<	real(kind=dp) :: listSkin = 1.0
664	>
665
666		!! initialize local variables
667
#	Line 574 \| Line 751 \| contains
751		iend = nGroups - 1
752		#endif
753		outer: do i = istart, iend
577	–
578	–	#ifdef IS_MPI
579	–	me_i = atid_row(i)
580	–	#else
581	–	me_i = atid(i)
582	–	#endif
754
755		if (update_nlist) point(i) = nlist + 1
756
#	Line 615 \| Line 786 \| contains
786		me_j = atid(j)
787		call get_interatomic_vector(q_group(:,i), &
788		q_group(:,j), d_grp, rgrpsq)
789	<	#endif
789	>	#endif
790
791	<	if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
791	>	if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
792		if (update_nlist) then
793		nlist = nlist + 1
794
#	Line 738 \| Line 909 \| contains
909		endif
910		end if
911		enddo
912	+
913		enddo outer
914
915		if (update_nlist) then
#	Line 817 \| Line 989 \| contains
989
990		if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
991
992	<	if (FF_uses_RF .and. SIM_uses_RF) then
992	>	if (electrostaticSummationMethod == REACTION_FIELD) then
993
994		#ifdef IS_MPI
995		call scatter(rf_Row,rf,plan_atom_row_3d)
#	Line 901 \| Line 1073 \| contains
1073		real ( kind = dp ), intent(inout) :: rijsq
1074		real ( kind = dp ) :: r
1075		real ( kind = dp ), intent(inout) :: d(3)
904	–	real ( kind = dp ) :: ebalance
1076		integer :: me_i, me_j
1077
1078		integer :: iHash
#	Line 928 \| Line 1099 \| contains
1099		call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1100		pot, eFrame, f, t, do_pot)
1101
1102	<	if (FF_uses_RF .and. SIM_uses_RF) then
1102	>	if (electrostaticSummationMethod == REACTION_FIELD) then
1103
1104		! CHECK ME (RF needs to know about all electrostatic types)
1105		call accumulate_rf(i, j, r, eFrame, sw)
#	Line 991 \| Line 1162 \| contains
1162
1163		integer :: me_i, me_j, iHash
1164
1165	+	r = sqrt(rijsq)
1166	+
1167		#ifdef IS_MPI
1168		me_i = atid_row(i)
1169		me_j = atid_col(j)
#	Line 1208 \| Line 1381 \| contains
1381
1382		function FF_RequiresPostpairCalc() result(doesit)
1383		logical :: doesit
1384	<	doesit = FF_uses_RF
1384	>	if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1385		end function FF_RequiresPostpairCalc
1386
1387		#ifdef PROFILE

Diff Legend

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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents): Revision 2270 by gezelter, Tue Aug 9 22:33:37 2005 UTC vs. Revision 2342 by chrisfen, Tue Oct 4 19:32:58 2005 UTC

Diff Legend

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2270 by gezelter, Tue Aug 9 22:33:37 2005 UTC vs.
Revision 2342 by chrisfen, Tue Oct 4 19:32:58 2005 UTC