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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2309 by chrisfen, Sun Sep 18 20:45:38 2005 UTC vs.
Revision 2375 by gezelter, Mon Oct 17 19:12:45 2005 UTC

#	Line 45 \| Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.46 2005-09-18 20:45:38 chrisfen Exp $, $Date: 2005-09-18 20:45:38 $, $Name: not supported by cvs2svn $, $Revision: 1.46 $
48	>	!! @version $Id: doForces.F90,v 1.59 2005-10-17 19:12:34 gezelter Exp $, $Date: 2005-10-17 19:12:34 $, $Name: not supported by cvs2svn $, $Revision: 1.59 $
49
50
51		module doForces
#	Line 59 \| Line 59 \| module doForces
59		use sticky
60		use electrostatic_module
61		use reaction_field_module
62	<	use gb_pair
62	>	use gayberne
63		use shapes
64		use vector_class
65		use eam
#	Line 75 \| Line 75 \| module doForces
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
#	Line 85 \| 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
#	Line 122 \| Line 124 \| module doForces
124		! Bit hash to determine pair-pair interactions.
125		integer, dimension(:,:), allocatable :: InteractionHash
126		real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
127	<	real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
128	<	integer, dimension(:), allocatable :: groupToGtype
129	<	real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
127	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
128	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
129	>
130	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
131	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
132	>
133	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
134	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
135		type ::gtypeCutoffs
136		real(kind=dp) :: rcut
137		real(kind=dp) :: rcutsq
#	Line 134 \| Line 141 \| module doForces
141
142		integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
143		real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
144	+	real(kind=dp),save :: listSkin
145
146		contains
147
#	Line 177 \| Line 185 \| contains
185
186		if (.not. allocated(InteractionHash)) then
187		allocate(InteractionHash(nAtypes,nAtypes))
188	+	else
189	+	deallocate(InteractionHash)
190	+	allocate(InteractionHash(nAtypes,nAtypes))
191		endif
192
193		if (.not. allocated(atypeMaxCutoff)) then
194		allocate(atypeMaxCutoff(nAtypes))
195	+	else
196	+	deallocate(atypeMaxCutoff)
197	+	allocate(atypeMaxCutoff(nAtypes))
198		endif
199
200		do i = 1, nAtypes
#	Line 257 \| Line 271 \| contains
271		logical :: GtypeFound
272
273		integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend
274	<	integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
274	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
275		integer :: nGroupsInRow
276	<	real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
276	>	integer :: nGroupsInCol
277	>	integer :: nGroupTypesRow,nGroupTypesCol
278	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
279		real(kind=dp) :: biggestAtypeCutoff
280
281		stat = 0
#	Line 273 \| Line 289 \| contains
289		endif
290		#ifdef IS_MPI
291		nGroupsInRow = getNgroupsInRow(plan_group_row)
292	+	nGroupsInCol = getNgroupsInCol(plan_group_col)
293		#endif
294		nAtypes = getSize(atypes)
295		! Set all of the initial cutoffs to zero.
#	Line 288 \| Line 305 \| contains
305		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
306
307
308	<	if (i_is_LJ) then
309	<	thisRcut = getSigma(i) * 2.5_dp
310	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
311	<	endif
312	<	if (i_is_Elect) then
313	<	thisRcut = defaultRcut
314	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
315	<	endif
316	<	if (i_is_Sticky) then
317	<	thisRcut = getStickyCut(i)
318	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
319	<	endif
320	<	if (i_is_StickyP) then
321	<	thisRcut = getStickyPowerCut(i)
322	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308	>	if (haveDefaultCutoffs) then
309	>	atypeMaxCutoff(i) = defaultRcut
310	>	else
311	>	if (i_is_LJ) then
312	>	thisRcut = getSigma(i) * 2.5_dp
313	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314	>	endif
315	>	if (i_is_Elect) then
316	>	thisRcut = defaultRcut
317	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
318	>	endif
319	>	if (i_is_Sticky) then
320	>	thisRcut = getStickyCut(i)
321	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322	>	endif
323	>	if (i_is_StickyP) then
324	>	thisRcut = getStickyPowerCut(i)
325	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326	>	endif
327	>	if (i_is_GB) then
328	>	thisRcut = getGayBerneCut(i)
329	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330	>	endif
331	>	if (i_is_EAM) then
332	>	thisRcut = getEAMCut(i)
333	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334	>	endif
335	>	if (i_is_Shape) then
336	>	thisRcut = getShapeCut(i)
337	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
338	>	endif
339		endif
307	–	if (i_is_GB) then
308	–	thisRcut = getGayBerneCut(i)
309	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
310	–	endif
311	–	if (i_is_EAM) then
312	–	thisRcut = getEAMCut(i)
313	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314	–	endif
315	–	if (i_is_Shape) then
316	–	thisRcut = getShapeCut(i)
317	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
318	–	endif
340
341	+
342		if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
343		biggestAtypeCutoff = atypeMaxCutoff(i)
344		endif
345	+
346		endif
347		enddo
348
349	<	nGroupTypes = 0
349	>
350
351		istart = 1
352	+	jstart = 1
353		#ifdef IS_MPI
354		iend = nGroupsInRow
355	+	jend = nGroupsInCol
356		#else
357		iend = nGroups
358	+	jend = nGroups
359		#endif
360
361		!! allocate the groupToGtype and gtypeMaxCutoff here.
362	<	if(.not.allocated(groupToGtype)) then
363	<	allocate(groupToGtype(iend))
364	<	allocate(groupMaxCutoff(iend))
365	<	allocate(gtypeMaxCutoff(iend))
366	<	groupMaxCutoff = 0.0_dp
367	<	gtypeMaxCutoff = 0.0_dp
362	>	if(.not.allocated(groupToGtypeRow)) then
363	>	! allocate(groupToGtype(iend))
364	>	allocate(groupToGtypeRow(iend))
365	>	else
366	>	deallocate(groupToGtypeRow)
367	>	allocate(groupToGtypeRow(iend))
368		endif
369	+	if(.not.allocated(groupMaxCutoffRow)) then
370	+	allocate(groupMaxCutoffRow(iend))
371	+	else
372	+	deallocate(groupMaxCutoffRow)
373	+	allocate(groupMaxCutoffRow(iend))
374	+	end if
375	+
376	+	if(.not.allocated(gtypeMaxCutoffRow)) then
377	+	allocate(gtypeMaxCutoffRow(iend))
378	+	else
379	+	deallocate(gtypeMaxCutoffRow)
380	+	allocate(gtypeMaxCutoffRow(iend))
381	+	endif
382	+
383	+
384	+	#ifdef IS_MPI
385	+	! We only allocate new storage if we are in MPI because Ncol /= Nrow
386	+	if(.not.associated(groupToGtypeCol)) then
387	+	allocate(groupToGtypeCol(jend))
388	+	else
389	+	deallocate(groupToGtypeCol)
390	+	allocate(groupToGtypeCol(jend))
391	+	end if
392	+
393	+	if(.not.associated(groupToGtypeCol)) then
394	+	allocate(groupToGtypeCol(jend))
395	+	else
396	+	deallocate(groupToGtypeCol)
397	+	allocate(groupToGtypeCol(jend))
398	+	end if
399	+	if(.not.associated(gtypeMaxCutoffCol)) then
400	+	allocate(gtypeMaxCutoffCol(jend))
401	+	else
402	+	deallocate(gtypeMaxCutoffCol)
403	+	allocate(gtypeMaxCutoffCol(jend))
404	+	end if
405	+
406	+	groupMaxCutoffCol = 0.0_dp
407	+	gtypeMaxCutoffCol = 0.0_dp
408	+
409	+	#endif
410	+	groupMaxCutoffRow = 0.0_dp
411	+	gtypeMaxCutoffRow = 0.0_dp
412	+
413	+
414		!! first we do a single loop over the cutoff groups to find the
415		!! largest cutoff for any atypes present in this group. We also
416		!! create gtypes at this point.
417
418		tol = 1.0d-6
419	<
419	>	nGroupTypesRow = 0
420	>
421		do i = istart, iend
422		n_in_i = groupStartRow(i+1) - groupStartRow(i)
423	<	groupMaxCutoff(i) = 0.0_dp
423	>	groupMaxCutoffRow(i) = 0.0_dp
424		do ia = groupStartRow(i), groupStartRow(i+1)-1
425		atom1 = groupListRow(ia)
426		#ifdef IS_MPI
#	Line 356 \| Line 428 \| contains
428		#else
429		me_i = atid(atom1)
430		#endif
431	<	if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
432	<	groupMaxCutoff(i)=atypeMaxCutoff(me_i)
431	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
432	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
433		endif
434		enddo
435
436	<	if (nGroupTypes.eq.0) then
437	<	nGroupTypes = nGroupTypes + 1
438	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
439	<	groupToGtype(i) = nGroupTypes
436	>	if (nGroupTypesRow.eq.0) then
437	>	nGroupTypesRow = nGroupTypesRow + 1
438	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
439	>	groupToGtypeRow(i) = nGroupTypesRow
440		else
441		GtypeFound = .false.
442	<	do g = 1, nGroupTypes
443	<	if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
444	<	groupToGtype(i) = g
442	>	do g = 1, nGroupTypesRow
443	>	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
444	>	groupToGtypeRow(i) = g
445		GtypeFound = .true.
446		endif
447		enddo
448		if (.not.GtypeFound) then
449	<	nGroupTypes = nGroupTypes + 1
450	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
451	<	groupToGtype(i) = nGroupTypes
449	>	nGroupTypesRow = nGroupTypesRow + 1
450	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
451	>	groupToGtypeRow(i) = nGroupTypesRow
452		endif
453		endif
454		enddo
455
456	+	#ifdef IS_MPI
457	+	do j = jstart, jend
458	+	n_in_j = groupStartCol(j+1) - groupStartCol(j)
459	+	groupMaxCutoffCol(j) = 0.0_dp
460	+	do ja = groupStartCol(j), groupStartCol(j+1)-1
461	+	atom1 = groupListCol(ja)
462	+
463	+	me_j = atid_col(atom1)
464	+
465	+	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
466	+	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
467	+	endif
468	+	enddo
469	+
470	+	if (nGroupTypesCol.eq.0) then
471	+	nGroupTypesCol = nGroupTypesCol + 1
472	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
473	+	groupToGtypeCol(j) = nGroupTypesCol
474	+	else
475	+	GtypeFound = .false.
476	+	do g = 1, nGroupTypesCol
477	+	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
478	+	groupToGtypeCol(j) = g
479	+	GtypeFound = .true.
480	+	endif
481	+	enddo
482	+	if (.not.GtypeFound) then
483	+	nGroupTypesCol = nGroupTypesCol + 1
484	+	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
485	+	groupToGtypeCol(j) = nGroupTypesCol
486	+	endif
487	+	endif
488	+	enddo
489	+
490	+	#else
491	+	! Set pointers to information we just found
492	+	nGroupTypesCol = nGroupTypesRow
493	+	groupToGtypeCol => groupToGtypeRow
494	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
495	+	groupMaxCutoffCol => groupMaxCutoffRow
496	+	#endif
497	+
498	+
499	+
500	+
501	+
502		!! allocate the gtypeCutoffMap here.
503	<	allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
503	>	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
504		!! then we do a double loop over all the group TYPES to find the cutoff
505		!! map between groups of two types
506	<
507	<	do i = 1, nGroupTypes
508	<	do j = 1, nGroupTypes
506	>	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507	>
508	>	do i = 1, nGroupTypesRow
509	>	do j = 1, nGroupTypesCol
510
511		select case(cutoffPolicy)
512		case(TRADITIONAL_CUTOFF_POLICY)
513	<	thisRcut = maxval(gtypeMaxCutoff)
513	>	thisRcut = tradRcut
514		case(MIX_CUTOFF_POLICY)
515	<	thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
515	>	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516		case(MAX_CUTOFF_POLICY)
517	<	thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
517	>	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
518		case default
519		call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520		return
#	Line 403 \| Line 522 \| contains
522		gtypeCutoffMap(i,j)%rcut = thisRcut
523		gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
524		skin = defaultRlist - defaultRcut
525	+	listSkin = skin ! set neighbor list skin thickness
526		gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
527
528	+	! sanity check
529	+
530	+	if (haveDefaultCutoffs) then
531	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
532	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
533	+	endif
534	+	endif
535		enddo
536		enddo
537	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
538	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
539	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
540	+	#ifdef IS_MPI
541	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
542	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
543	+	#endif
544	+	groupMaxCutoffCol => null()
545	+	gtypeMaxCutoffCol => null()
546
547		haveGtypeCutoffMap = .true.
548		end subroutine createGtypeCutoffMap
#	Line 419 \| Line 555 \| contains
555		defaultRsw = defRsw
556		defaultRlist = defRlist
557		cutoffPolicy = cutPolicy
558	+
559	+	haveDefaultCutoffs = .true.
560		end subroutine setDefaultCutoffs
561
562		subroutine setCutoffPolicy(cutPolicy)
#	Line 504 \| Line 642 \| contains
642		subroutine init_FF(thisESM, thisStat)
643
644		integer, intent(in) :: thisESM
507	–	real(kind=dp), intent(in) :: dampingAlpha
645		integer, intent(out) :: thisStat
646		integer :: my_status, nMatches
647		integer, pointer :: MatchList(:) => null()
#	Line 547 \| Line 684 \| contains
684		!! check to make sure the reaction field setting makes sense
685
686		if (FF_uses_Dipoles) then
687	<	if (electrostaticSummationMethod == 3) then
687	>	if (electrostaticSummationMethod == REACTION_FIELD) then
688		dielect = getDielect()
689		call initialize_rf(dielect)
690		endif
691		else
692	<	if (electrostaticSummationMethod == 3) then
692	>	if (electrostaticSummationMethod == REACTION_FIELD) then
693		write(default_error,*) 'Using Reaction Field with no dipoles? Huh?'
694		thisStat = -1
695		haveSaneForceField = .false.
#	Line 568 \| Line 705 \| contains
705		haveSaneForceField = .false.
706		return
707		end if
571	–	endif
572	–
573	–	if (FF_uses_GayBerne) then
574	–	call check_gb_pair_FF(my_status)
575	–	if (my_status .ne. 0) then
576	–	thisStat = -1
577	–	haveSaneForceField = .false.
578	–	return
579	–	endif
708		endif
709
710		if (.not. haveNeighborList) then
#	Line 612 \| Line 740 \| contains
740
741		!! Stress Tensor
742		real( kind = dp), dimension(9) :: tau
743	<	real ( kind = dp ) :: pot
743	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
744		logical ( kind = 2) :: do_pot_c, do_stress_c
745		logical :: do_pot
746		logical :: do_stress
747		logical :: in_switching_region
748		#ifdef IS_MPI
749	<	real( kind = DP ) :: pot_local
749	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
750		integer :: nAtomsInRow
751		integer :: nAtomsInCol
752		integer :: nprocs
#	Line 643 \| Line 771 \| contains
771		integer :: propPack_i, propPack_j
772		integer :: loopStart, loopEnd, loop
773		integer :: iHash
774	<	real(kind=dp) :: listSkin = 1.0
774	>
775
776		!! initialize local variables
777
#	Line 768 \| Line 896 \| contains
896		me_j = atid(j)
897		call get_interatomic_vector(q_group(:,i), &
898		q_group(:,j), d_grp, rgrpsq)
899	<	#endif
899	>	#endif
900
901	<	if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
901	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
902		if (update_nlist) then
903		nlist = nlist + 1
904
#	Line 891 \| Line 1019 \| contains
1019		endif
1020		end if
1021		enddo
1022	+
1023		enddo outer
1024
1025		if (update_nlist) then
#	Line 950 \| Line 1079 \| contains
1079
1080		if (do_pot) then
1081		! scatter/gather pot_row into the members of my column
1082	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1083	<
1082	>	do i = 1,LR_POT_TYPES
1083	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1084	>	end do
1085		! scatter/gather pot_local into all other procs
1086		! add resultant to get total pot
1087		do i = 1, nlocal
1088	<	pot_local = pot_local + pot_Temp(i)
1088	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1089	>	+ pot_Temp(1:LR_POT_TYPES,i)
1090		enddo
1091
1092		pot_Temp = 0.0_DP
1093	<
1094	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1093	>	do i = 1,LR_POT_TYPES
1094	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1095	>	end do
1096		do i = 1, nlocal
1097	<	pot_local = pot_local + pot_Temp(i)
1097	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1098	>	+ pot_Temp(1:LR_POT_TYPES,i)
1099		enddo
1100
1101		endif
#	Line 970 \| Line 1103 \| contains
1103
1104		if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1105
1106	<	if (electrostaticSummationMethod == 3) then
1106	>	if (electrostaticSummationMethod == REACTION_FIELD) then
1107
1108		#ifdef IS_MPI
1109		call scatter(rf_Row,rf,plan_atom_row_3d)
#	Line 1001 \| Line 1134 \| contains
1134		!! potential and torques:
1135		call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1136		#ifdef IS_MPI
1137	<	pot_local = pot_local + rfpot
1137	>	pot_local(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
1138		#else
1139	<	pot = pot + rfpot
1139	>	pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1140
1141		#endif
1142		endif
#	Line 1015 \| Line 1148 \| contains
1148		#ifdef IS_MPI
1149
1150		if (do_pot) then
1151	<	pot = pot + pot_local
1151	>	pot(1:LR_POT_TYPES) = pot(1:LR_POT_TYPES) &
1152	>	+ pot_local(1:LR_POT_TYPES)
1153		!! we assume the c code will do the allreduce to get the total potential
1154		!! we could do it right here if we needed to...
1155		endif
#	Line 1041 \| Line 1175 \| contains
1175		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1176		eFrame, A, f, t, pot, vpair, fpair)
1177
1178	<	real( kind = dp ) :: pot, vpair, sw
1178	>	real( kind = dp ) :: vpair, sw
1179	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1180		real( kind = dp ), dimension(3) :: fpair
1181		real( kind = dp ), dimension(nLocal) :: mfact
1182		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1073 \| Line 1208 \| contains
1208		iHash = InteractionHash(me_i, me_j)
1209
1210		if ( iand(iHash, LJ_PAIR).ne.0 ) then
1211	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1211	>	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1212	>	pot(VDW_POT), f, do_pot)
1213		endif
1214
1215		if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1216		call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1217	<	pot, eFrame, f, t, do_pot)
1217	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1218
1219	<	if (electrostaticSummationMethod == 3) then
1219	>	if (electrostaticSummationMethod == REACTION_FIELD) then
1220
1221		! CHECK ME (RF needs to know about all electrostatic types)
1222		call accumulate_rf(i, j, r, eFrame, sw)
#	Line 1091 \| Line 1227 \| contains
1227
1228		if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1229		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1230	<	pot, A, f, t, do_pot)
1230	>	pot(HB_POT), A, f, t, do_pot)
1231		endif
1232
1233		if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1234		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1235	<	pot, A, f, t, do_pot)
1235	>	pot(HB_POT), A, f, t, do_pot)
1236		endif
1237
1238		if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1239		call do_gb_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(iHash, GAYBERNE_LJ).ne.0 ) then
1244	<	! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1245	<	! pot, A, f, t, do_pot)
1244	>	call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1245	>	pot(VDW_POT), A, f, t, do_pot)
1246		endif
1247
1248		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1249	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1250	<	do_pot)
1249	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1250	>	pot(METALLIC_POT), f, do_pot)
1251		endif
1252
1253		if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1254		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1255	<	pot, A, f, t, do_pot)
1255	>	pot(VDW_POT), A, f, t, do_pot)
1256		endif
1257
1258		if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1259		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1260	<	pot, A, f, t, do_pot)
1260	>	pot(VDW_POT), A, f, t, do_pot)
1261		endif
1262
1263		end subroutine do_pair
#	Line 1129 \| Line 1265 \| contains
1265		subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1266		do_pot, do_stress, eFrame, A, f, t, pot)
1267
1268	<	real( kind = dp ) :: pot, sw
1268	>	real( kind = dp ) :: sw
1269	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1270		real( kind = dp ), dimension(9,nLocal) :: eFrame
1271		real (kind=dp), dimension(9,nLocal) :: A
1272		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1164 \| Line 1301 \| contains
1301
1302		subroutine do_preforce(nlocal,pot)
1303		integer :: nlocal
1304	<	real( kind = dp ) :: pot
1304	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1305
1306		if (FF_uses_EAM .and. SIM_uses_EAM) then
1307	<	call calc_EAM_preforce_Frho(nlocal,pot)
1307	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1308		endif
1309
1310
#	Line 1362 \| Line 1499 \| contains
1499
1500		function FF_RequiresPostpairCalc() result(doesit)
1501		logical :: doesit
1502	<	if (electrostaticSummationMethod == 3) doesit = .true.
1502	>	if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1503		end function FF_RequiresPostpairCalc
1504
1505		#ifdef PROFILE

Diff Legend

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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents): Revision 2309 by chrisfen, Sun Sep 18 20:45:38 2005 UTC vs. Revision 2375 by gezelter, Mon Oct 17 19:12:45 2005 UTC

Diff Legend

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2309 by chrisfen, Sun Sep 18 20:45:38 2005 UTC vs.
Revision 2375 by gezelter, Mon Oct 17 19:12:45 2005 UTC