[ViewVC] Diff of: group/trunk/OOPSE/libmdtools/do

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 898 by chuckv, Mon Jan 5 22:49:14 2004 UTC vs.
Revision 941 by gezelter, Tue Jan 13 23:01:43 2004 UTC

#	Line 4 \| Line 4
4
5		!! @author Charles F. Vardeman II
6		!! @author Matthew Meineke
7	<	!! @version $Id: do_Forces.F90,v 1.44 2004-01-05 22:49:14 chuckv Exp $, $Date: 2004-01-05 22:49:14 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $
7	>	!! @version $Id: do_Forces.F90,v 1.46 2004-01-13 23:01:43 gezelter Exp $, $Date: 2004-01-13 23:01:43 $, $Name: not supported by cvs2svn $, $Revision: 1.46 $
8
9		module do_Forces
10		use force_globals
#	Line 15 \| Line 15 \| module do_Forces
15		use lj
16		use sticky_pair
17		use dipole_dipole
18	+	use charge_charge
19		use reaction_field
20		use gb_pair
21		use vector_class
#	Line 30 \| Line 31 \| module do_Forces
31		#define __FORTRAN90
32		#include "fForceField.h"
33
34	<	logical, save :: do_forces_initialized = .false., haveRlist = .false.
34	>	logical, save :: haveRlist = .false.
35	>	logical, save :: haveNeighborList = .false.
36		logical, save :: havePolicies = .false.
37	+	logical, save :: haveSIMvariables = .false.
38	+	logical, save :: havePropertyMap = .false.
39	+	logical, save :: haveSaneForceField = .false.
40		logical, save :: FF_uses_LJ
41		logical, save :: FF_uses_sticky
42	+	logical, save :: FF_uses_charges
43		logical, save :: FF_uses_dipoles
44		logical, save :: FF_uses_RF
45		logical, save :: FF_uses_GB
46		logical, save :: FF_uses_EAM
47	+	logical, save :: SIM_uses_LJ
48	+	logical, save :: SIM_uses_sticky
49	+	logical, save :: SIM_uses_charges
50	+	logical, save :: SIM_uses_dipoles
51	+	logical, save :: SIM_uses_RF
52	+	logical, save :: SIM_uses_GB
53	+	logical, save :: SIM_uses_EAM
54	+	logical, save :: SIM_requires_postpair_calc
55	+	logical, save :: SIM_requires_prepair_calc
56	+	logical, save :: SIM_uses_directional_atoms
57	+	logical, save :: SIM_uses_PBC
58
59		real(kind=dp), save :: rlist, rlistsq
60
#	Line 45 \| Line 62 \| module do_Forces
62		public :: do_force_loop
63		public :: setRlistDF
64
65	+
66		#ifdef PROFILE
67		public :: getforcetime
68		real, save :: forceTime = 0
#	Line 52 \| Line 70 \| module do_Forces
70		integer :: nLoops
71		#endif
72
73	<	logical, allocatable :: propertyMapI(:,:)
74	<	logical, allocatable :: propertyMapJ(:,:)
73	>	type :: Properties
74	>	logical :: is_lj = .false.
75	>	logical :: is_sticky = .false.
76	>	logical :: is_dp = .false.
77	>	logical :: is_gb = .false.
78	>	logical :: is_eam = .false.
79	>	logical :: is_charge = .false.
80	>	real(kind=DP) :: charge = 0.0_DP
81	>	real(kind=DP) :: dipole_moment = 0.0_DP
82	>	end type Properties
83
84	+	type(Properties), dimension(:),allocatable :: PropertyMap
85	+
86		contains
87
88		subroutine setRlistDF( this_rlist )
#	Line 65 \| Line 93 \| contains
93		rlistsq = rlist * rlist
94
95		haveRlist = .true.
68	–	if( havePolicies ) do_forces_initialized = .true.
96
97		end subroutine setRlistDF
98
99	+	subroutine createPropertyMap(status)
100	+	integer :: nAtypes
101	+	integer :: status
102	+	integer :: i
103	+	logical :: thisProperty
104	+	real (kind=DP) :: thisDPproperty
105	+
106	+	status = 0
107	+
108	+	nAtypes = getSize(atypes)
109	+
110	+	if (nAtypes == 0) then
111	+	status = -1
112	+	return
113	+	end if
114	+
115	+	if (.not. allocated(PropertyMap)) then
116	+	allocate(PropertyMap(nAtypes))
117	+	endif
118	+
119	+	do i = 1, nAtypes
120	+	call getElementProperty(atypes, i, "is_LJ", thisProperty)
121	+	PropertyMap(i)%is_LJ = thisProperty
122	+
123	+	call getElementProperty(atypes, i, "is_Charge", thisProperty)
124	+	PropertyMap(i)%is_Charge = thisProperty
125	+
126	+	if (thisProperty) then
127	+	call getElementProperty(atypes, i, "charge", thisDPproperty)
128	+	PropertyMap(i)%charge = thisDPproperty
129	+	endif
130	+
131	+	call getElementProperty(atypes, i, "is_DP", thisProperty)
132	+	PropertyMap(i)%is_DP = thisProperty
133	+
134	+	if (thisProperty) then
135	+	call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
136	+	PropertyMap(i)%dipole_moment = thisDPproperty
137	+	endif
138	+
139	+	call getElementProperty(atypes, i, "is_Sticky", thisProperty)
140	+	PropertyMap(i)%is_Sticky = thisProperty
141	+	call getElementProperty(atypes, i, "is_GB", thisProperty)
142	+	PropertyMap(i)%is_GB = thisProperty
143	+	call getElementProperty(atypes, i, "is_EAM", thisProperty)
144	+	PropertyMap(i)%is_EAM = thisProperty
145	+	end do
146	+
147	+	havePropertyMap = .true.
148	+
149	+	end subroutine createPropertyMap
150	+
151	+	subroutine setSimVariables()
152	+	SIM_uses_LJ = SimUsesLJ()
153	+	SIM_uses_sticky = SimUsesSticky()
154	+	SIM_uses_charges = SimUsesCharges()
155	+	SIM_uses_dipoles = SimUsesDipoles()
156	+	SIM_uses_RF = SimUsesRF()
157	+	SIM_uses_GB = SimUsesGB()
158	+	SIM_uses_EAM = SimUsesEAM()
159	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
160	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
161	+	SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
162	+	SIM_uses_PBC = SimUsesPBC()
163	+
164	+	haveSIMvariables = .true.
165	+
166	+	return
167	+	end subroutine setSimVariables
168	+
169	+	subroutine doReadyCheck(error)
170	+	integer, intent(out) :: error
171	+
172	+	integer :: myStatus
173	+
174	+	error = 0
175	+
176	+	if (.not. havePropertyMap) then
177	+
178	+	myStatus = 0
179	+
180	+	call createPropertyMap(myStatus)
181	+
182	+	if (myStatus .ne. 0) then
183	+	write(default_error, *) 'createPropertyMap failed in do_Forces!'
184	+	error = -1
185	+	return
186	+	endif
187	+	endif
188	+
189	+	if (.not. haveSIMvariables) then
190	+	call setSimVariables()
191	+	endif
192	+
193	+	if (.not. haveRlist) then
194	+	write(default_error, *) 'rList has not been set in do_Forces!'
195	+	error = -1
196	+	return
197	+	endif
198	+
199	+	if (SIM_uses_LJ .and. FF_uses_LJ) then
200	+	if (.not. havePolicies) then
201	+	write(default_error, *) 'LJ mixing Policies have not been set in do_Forces!'
202	+	error = -1
203	+	return
204	+	endif
205	+	endif
206	+
207	+	if (.not. haveNeighborList) then
208	+	write(default_error, *) 'neighbor list has not been initialized in do_Forces!'
209	+	error = -1
210	+	return
211	+	end if
212	+
213	+	if (.not. haveSaneForceField) then
214	+	write(default_error, *) 'Force Field is not sane in do_Forces!'
215	+	error = -1
216	+	return
217	+	end if
218	+
219	+	#ifdef IS_MPI
220	+	if (.not. isMPISimSet()) then
221	+	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
222	+	error = -1
223	+	return
224	+	endif
225	+	#endif
226	+	return
227	+	end subroutine doReadyCheck
228	+
229	+
230		subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
231
232		integer, intent(in) :: LJMIXPOLICY
#	Line 93 \| Line 251 \| contains
251
252		FF_uses_LJ = .false.
253		FF_uses_sticky = .false.
254	+	FF_uses_charges = .false.
255		FF_uses_dipoles = .false.
256		FF_uses_GB = .false.
257		FF_uses_EAM = .false.
258
259		call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
260		if (nMatches .gt. 0) FF_uses_LJ = .true.
261	+
262	+	call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
263	+	if (nMatches .gt. 0) FF_uses_charges = .true.
264
265		call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
266		if (nMatches .gt. 0) FF_uses_dipoles = .true.
#	Line 113 \| Line 275 \| contains
275		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
276		if (nMatches .gt. 0) FF_uses_EAM = .true.
277
278	+	!! Assume sanity (for the sake of argument)
279	+	haveSaneForceField = .true.
280	+
281		!! check to make sure the FF_uses_RF setting makes sense
282
283		if (FF_uses_dipoles) then
#	Line 124 \| Line 289 \| contains
289		if (FF_uses_RF) then
290		write(default_error,*) 'Using Reaction Field with no dipoles? Huh?'
291		thisStat = -1
292	+	haveSaneForceField = .false.
293		return
294		endif
295		endif
#	Line 138 \| Line 304 \| contains
304		case default
305		write(default_error,*) 'unknown LJ Mixing Policy!'
306		thisStat = -1
307	+	haveSaneForceField = .false.
308		return
309		end select
310		if (my_status /= 0) then
311		thisStat = -1
312	+	haveSaneForceField = .false.
313		return
314		end if
315	+	havePolicies = .true.
316		endif
317
318		if (FF_uses_sticky) then
319		call check_sticky_FF(my_status)
320		if (my_status /= 0) then
321		thisStat = -1
322	+	haveSaneForceField = .false.
323		return
324		end if
325		endif
#	Line 158 \| Line 328 \| contains
328		if (FF_uses_EAM) then
329		call init_EAM_FF(my_status)
330		if (my_status /= 0) then
331	<	write(,) "init_EAM_FF returned a bad status"
331	>	write(default_error, *) "init_EAM_FF returned a bad status"
332		thisStat = -1
333	+	haveSaneForceField = .false.
334		return
335		end if
336		endif
337
167	–
168	–
338		if (FF_uses_GB) then
339		call check_gb_pair_FF(my_status)
340		if (my_status .ne. 0) then
341		thisStat = -1
342	+	haveSaneForceField = .false.
343		return
344		endif
345		endif
346
347		if (FF_uses_GB .and. FF_uses_LJ) then
348		endif
349	<	if (.not. do_forces_initialized) then
349	>	if (.not. haveNeighborList) then
350		!! Create neighbor lists
351		call expandNeighborList(nLocal, my_status)
352		if (my_Status /= 0) then
#	Line 184 \| Line 354 \| contains
354		thisStat = -1
355		return
356		endif
357	+	haveNeighborList = .true.
358		endif
359
189	–
190	–	havePolicies = .true.
191	–	if( haveRlist ) do_forces_initialized = .true.
192	–
360		end subroutine init_FF
361
362
#	Line 246 \| Line 413 \| contains
413		natoms = nlocal
414		#endif
415
416	<	call check_initialization(localError)
416	>	call doReadyCheck(localError)
417		if ( localError .ne. 0 ) then
418	<	call handleError("do_force_loop","Not Initialized")
418	>	call handleError("do_force_loop", "Not Initialized")
419		error = -1
420		return
421		end if
#	Line 257 \| Line 424 \| contains
424		do_pot = do_pot_c
425		do_stress = do_stress_c
426
260	–
261	–	#ifdef IS_MPI
262	–	if (.not.allocated(propertyMapI)) then
263	–	allocate(propertyMapI(5,nrow))
264	–	endif
265	–
266	–	do i = 1, nrow
267	–	me_i = atid_row(i)
268	–	#else
269	–	if (.not.allocated(propertyMapI)) then
270	–	allocate(propertyMapI(5,nlocal))
271	–	endif
272	–
273	–	do i = 1, natoms
274	–	me_i = atid(i)
275	–	#endif
276	–
277	–	propertyMapI(1:5,i) = .false.
278	–
279	–	call getElementProperty(atypes, me_i, "propertyPack", propPack_i)
280	–
281	–	! unpack the properties
282	–
283	–	if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
284	–	propertyMapI(1, i) = .true.
285	–	if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
286	–	propertyMapI(2, i) = .true.
287	–	if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
288	–	propertyMapI(3, i) = .true.
289	–	if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
290	–	propertyMapI(4, i) = .true.
291	–	if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
292	–	propertyMapI(5, i) = .true.
293	–
294	–	end do
295	–
296	–	#ifdef IS_MPI
297	–	if (.not.allocated(propertyMapJ)) then
298	–	allocate(propertyMapJ(5,ncol))
299	–	endif
300	–
301	–	do j = 1, ncol
302	–	me_j = atid_col(j)
303	–	#else
304	–	if (.not.allocated(propertyMapJ)) then
305	–	allocate(propertyMapJ(5,nlocal))
306	–	endif
307	–
308	–	do j = 1, natoms
309	–	me_j = atid(j)
310	–	#endif
311	–
312	–	propertyMapJ(1:5,j) = .false.
313	–
314	–	call getElementProperty(atypes, me_j, "propertyPack", propPack_j)
315	–
316	–	! unpack the properties
317	–
318	–	if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
319	–	propertyMapJ(1, j) = .true.
320	–	if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
321	–	propertyMapJ(2, j) = .true.
322	–	if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
323	–	propertyMapJ(3, j) = .true.
324	–	if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
325	–	propertyMapJ(4, j) = .true.
326	–	if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
327	–	propertyMapJ(5, j) = .true.
328	–
329	–	end do
330	–
427		! Gather all information needed by all force loops:
428
429		#ifdef IS_MPI
#	Line 335 \| Line 431 \| contains
431		call gather(q,q_Row,plan_row3d)
432		call gather(q,q_Col,plan_col3d)
433
434	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
434	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
435		call gather(u_l,u_l_Row,plan_row3d)
436		call gather(u_l,u_l_Col,plan_col3d)
437
#	Line 351 \| Line 447 \| contains
447		nloops = nloops + 1
448		#endif
449
450	<	if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
450	>	if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
451		!! See if we need to update neighbor lists
452		call checkNeighborList(nlocal, q, listSkin, update_nlist)
453		!! if_mpi_gather_stuff_for_prepair
#	Line 672 \| Line 768 \| contains
768		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
769		end do
770
771	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
771	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
772		t_temp = 0.0_dp
773		call scatter(t_Row,t_temp,plan_row3d)
774		do i = 1,nlocal
#	Line 706 \| Line 802 \| contains
802		endif
803		#endif
804
805	<	if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
805	>	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
806
807	<	if (FF_uses_RF .and. SimUsesRF()) then
807	>	if (FF_uses_RF .and. SIM_uses_RF) then
808
809		#ifdef IS_MPI
810		call scatter(rf_Row,rf,plan_row3d)
#	Line 726 \| Line 822 \| contains
822		#else
823		me_i = atid(i)
824		#endif
825	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
826	<	if ( is_DP_i ) then
827	<	call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
825	>
826	>	if (PropertyMap(me_i)%is_DP) then
827	>
828	>	mu_i = PropertyMap(me_i)%dipole_moment
829	>
830		!! The reaction field needs to include a self contribution
831		!! to the field:
832	<	call accumulate_self_rf(i, mu_i, u_l)
832	>	call accumulate_self_rf(i, mu_i, u_l)
833		!! Get the reaction field contribution to the
834		!! potential and torques:
835		call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
#	Line 813 \| Line 911 \| contains
911
912		#endif
913
914	<	if (FF_uses_LJ .and. SimUsesLJ()) then
914	>	if (FF_uses_LJ .and. SIM_uses_LJ) then
915
916	<	if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) &
916	>	if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) &
917		call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
918
919		endif
920
921	<	if (FF_uses_dipoles .and. SimUsesDipoles()) then
921	>	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
922
923	<	if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then
923	>	if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
924		call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
925		do_pot, do_stress)
926	<	if (FF_uses_RF .and. SimUsesRF()) then
926	>	if (FF_uses_RF .and. SIM_uses_RF) then
927		call accumulate_rf(i, j, r, u_l)
928		call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
929		endif
#	Line 833 \| Line 931 \| contains
931		endif
932		endif
933
934	<	if (FF_uses_Sticky .and. SimUsesSticky()) then
934	>	if (FF_uses_Sticky .and. SIM_uses_sticky) then
935
936	<	if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then
936	>	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
937		call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
938		do_pot, do_stress)
939		endif
940		endif
941
942
943	<	if (FF_uses_GB .and. SimUsesGB()) then
943	>	if (FF_uses_GB .and. SIM_uses_GB) then
944
945	<	if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then
945	>	if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
946		call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
947		do_pot, do_stress)
948		endif
#	Line 853 \| Line 951 \| contains
951
952
953
954	<	if (FF_uses_EAM .and. SimUsesEAM()) then
954	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
955
956	<	if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then
956	>	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
957		call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
958		endif
959
#	Line 900 \| Line 998 \| contains
998
999		#endif
1000
1001	<	if (FF_uses_EAM .and. SimUsesEAM()) then
1002	<	call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
1003	<	call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
906	<
907	<	if ( is_EAM_i .and. is_EAM_j ) &
1001	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1002	>
1003	>	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1004		call calc_EAM_prepair_rho(i, j, d, r, rijsq )
909	–	endif
1005
1006	+	endif
1007	+
1008		end subroutine do_prepair
1009
1010
#	Line 917 \| Line 1014 \| contains
1014		integer :: nlocal
1015		real( kind = dp ) :: pot
1016
1017	<	if (FF_uses_EAM .and. SimUsesEAM()) then
1017	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1018		call calc_EAM_preforce_Frho(nlocal,pot)
1019		endif
1020
#	Line 936 \| Line 1033 \| contains
1033		d(1:3) = q_j(1:3) - q_i(1:3)
1034
1035		! Wrap back into periodic box if necessary
1036	<	if ( SimUsesPBC() ) then
1036	>	if ( SIM_uses_PBC ) then
1037
1038		if( .not.boxIsOrthorhombic ) then
1039		! calc the scaled coordinates.
#	Line 976 \| Line 1073 \| contains
1073
1074		end subroutine get_interatomic_vector
1075
979	–	subroutine check_initialization(error)
980	–	integer, intent(out) :: error
981	–
982	–	error = 0
983	–	! Make sure we are properly initialized.
984	–	if (.not. do_forces_initialized) then
985	–	write(,) "Forces not initialized"
986	–	error = -1
987	–	return
988	–	endif
989	–
990	–	#ifdef IS_MPI
991	–	if (.not. isMPISimSet()) then
992	–	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
993	–	error = -1
994	–	return
995	–	endif
996	–	#endif
997	–
998	–	return
999	–	end subroutine check_initialization
1000	–
1001	–
1076		subroutine zero_work_arrays()
1077
1078		#ifdef IS_MPI
#	Line 1031 \| Line 1105 \| contains
1105		#endif
1106
1107
1108	<	if (FF_uses_EAM .and. SimUsesEAM()) then
1108	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1109		call clean_EAM()
1110		endif
1111

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents): Revision 898 by chuckv, Mon Jan 5 22:49:14 2004 UTC vs. Revision 941 by gezelter, Tue Jan 13 23:01:43 2004 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 898 by chuckv, Mon Jan 5 22:49:14 2004 UTC vs.
Revision 941 by gezelter, Tue Jan 13 23:01:43 2004 UTC