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

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 626 by mmeineke, Wed Jul 16 21:30:56 2003 UTC vs.
Revision 883 by chuckv, Thu Dec 18 20:46:45 2003 UTC

#	Line 4 \| Line 4
4
5		!! @author Charles F. Vardeman II
6		!! @author Matthew Meineke
7	<	!! @version $Id: do_Forces.F90,v 1.21 2003-07-16 21:30:55 mmeineke Exp $, $Date: 2003-07-16 21:30:55 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $
7	>	!! @version $Id: do_Forces.F90,v 1.38 2003-12-18 20:46:45 chuckv Exp $, $Date: 2003-12-18 20:46:45 $, $Name: not supported by cvs2svn $, $Revision: 1.38 $
8
9		module do_Forces
10		use force_globals
#	Line 18 \| Line 18 \| module do_Forces
18		use reaction_field
19		use gb_pair
20		use vector_class
21	+	use eam
22	+	use status
23		#ifdef IS_MPI
24		use mpiSimulation
25		#endif
#	Line 43 \| Line 45 \| contains
45		public :: do_force_loop
46		public :: setRlistDF
47
48	+	#ifdef PROFILE
49	+	public :: getforcetime
50	+	real, save :: forceTime = 0
51	+	real :: forceTimeInitial, forceTimeFinal
52	+	#endif
53	+
54		contains
55
56		subroutine setRlistDF( this_rlist )
#	Line 136 \| Line 144 \| contains
144
145		if (FF_uses_sticky) then
146		call check_sticky_FF(my_status)
147	+	if (my_status /= 0) then
148	+	thisStat = -1
149	+	return
150	+	end if
151	+	endif
152	+
153	+
154	+	if (FF_uses_EAM) then
155	+	call init_EAM_FF(my_status)
156		if (my_status /= 0) then
157	+	write(,) "init_EAM_FF returned a bad status"
158		thisStat = -1
159		return
160		end if
161		endif
162	+
163	+
164
165		if (FF_uses_GB) then
166		call check_gb_pair_FF(my_status)
#	Line 161 \| Line 181 \| contains
181		return
182		endif
183		endif
184	+
185
186		havePolicies = .true.
187		if( haveRlist ) do_forces_initialized = .true.
188	<
188	>
189		end subroutine init_FF
190
191
#	Line 192 \| Line 213 \| contains
213		real( kind = DP ) :: pot_local
214		integer :: nrow
215		integer :: ncol
216	+	integer :: nprocs
217		#endif
218		integer :: nlocal
219		integer :: natoms
#	Line 207 \| Line 229 \| contains
229		integer :: listerror, error
230		integer :: localError
231
232	<	real(kind=dp) :: listSkin = 1.0
211	<
232	>	real(kind=dp) :: listSkin = 1.0
233
234		!! initialize local variables
235
#	Line 221 \| Line 242 \| contains
242		nlocal = getNlocal()
243		natoms = nlocal
244		#endif
245	<
245	>
246		call check_initialization(localError)
247		if ( localError .ne. 0 ) then
248	+	call handleError("do_force_loop","Not Initialized")
249		error = -1
250		return
251		end if
#	Line 232 \| Line 254 \| contains
254		do_pot = do_pot_c
255		do_stress = do_stress_c
256
257	+
258		! Gather all information needed by all force loops:
259
260		#ifdef IS_MPI
#	Line 248 \| Line 271 \| contains
271		endif
272
273		#endif
274	<
274	>
275	>	!! Begin force loop timing:
276	>	#ifdef PROFILE
277	>	call cpu_time(forceTimeInitial)
278	>	nloops = nloops + 1
279	>	#endif
280	>
281		if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
282		!! See if we need to update neighbor lists
283		call checkNeighborList(nlocal, q, listSkin, update_nlist)
#	Line 256 \| Line 285 \| contains
285		!! do_prepair_loop_if_needed
286		!! if_mpi_scatter_stuff_from_prepair
287		!! if_mpi_gather_stuff_from_prepair_to_main_loop
288	+
289	+	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
290	+	#ifdef IS_MPI
291	+
292	+	if (update_nlist) then
293	+
294	+	!! save current configuration, construct neighbor list,
295	+	!! and calculate forces
296	+	call saveNeighborList(nlocal, q)
297	+
298	+	neighborListSize = size(list)
299	+	nlist = 0
300	+
301	+	do i = 1, nrow
302	+	point(i) = nlist + 1
303	+
304	+	prepair_inner: do j = 1, ncol
305	+
306	+	if (skipThisPair(i,j)) cycle prepair_inner
307	+
308	+	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
309	+
310	+	if (rijsq < rlistsq) then
311	+
312	+	nlist = nlist + 1
313	+
314	+	if (nlist > neighborListSize) then
315	+	call expandNeighborList(nlocal, listerror)
316	+	if (listerror /= 0) then
317	+	error = -1
318	+	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
319	+	return
320	+	end if
321	+	neighborListSize = size(list)
322	+	endif
323	+
324	+	list(nlist) = j
325	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)
326	+	endif
327	+	enddo prepair_inner
328	+	enddo
329	+
330	+	point(nrow + 1) = nlist + 1
331	+
332	+	else !! (of update_check)
333	+
334	+	! use the list to find the neighbors
335	+	do i = 1, nrow
336	+	JBEG = POINT(i)
337	+	JEND = POINT(i+1) - 1
338	+	! check thiat molecule i has neighbors
339	+	if (jbeg .le. jend) then
340	+
341	+	do jnab = jbeg, jend
342	+	j = list(jnab)
343	+
344	+	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
345	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
346	+	u_l, A, f, t, pot_local)
347	+
348	+	enddo
349	+	endif
350	+	enddo
351	+	endif
352	+
353	+	#else
354	+
355	+	if (update_nlist) then
356	+
357	+	! save current configuration, contruct neighbor list,
358	+	! and calculate forces
359	+	call saveNeighborList(natoms, q)
360	+
361	+	neighborListSize = size(list)
362	+
363	+	nlist = 0
364	+
365	+	do i = 1, natoms-1
366	+	point(i) = nlist + 1
367	+
368	+	prepair_inner: do j = i+1, natoms
369	+
370	+	if (skipThisPair(i,j)) cycle prepair_inner
371	+
372	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
373	+
374	+
375	+	if (rijsq < rlistsq) then
376	+
377	+
378	+	nlist = nlist + 1
379	+
380	+	if (nlist > neighborListSize) then
381	+	call expandNeighborList(natoms, listerror)
382	+	if (listerror /= 0) then
383	+	error = -1
384	+	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
385	+	return
386	+	end if
387	+	neighborListSize = size(list)
388	+	endif
389	+
390	+	list(nlist) = j
391	+
392	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
393	+	u_l, A, f, t, pot)
394	+
395	+	endif
396	+	enddo prepair_inner
397	+	enddo
398	+
399	+	point(natoms) = nlist + 1
400	+
401	+	else !! (update)
402	+
403	+	! use the list to find the neighbors
404	+	do i = 1, natoms-1
405	+	JBEG = POINT(i)
406	+	JEND = POINT(i+1) - 1
407	+	! check thiat molecule i has neighbors
408	+	if (jbeg .le. jend) then
409	+
410	+	do jnab = jbeg, jend
411	+	j = list(jnab)
412	+
413	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
414	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
415	+	u_l, A, f, t, pot)
416	+
417	+	enddo
418	+	endif
419	+	enddo
420	+	endif
421	+	#endif
422	+	!! Do rest of preforce calculations
423	+	!! do necessary preforce calculations
424	+	call do_preforce(nlocal,pot)
425	+	! we have already updated the neighbor list set it to false...
426	+	update_nlist = .false.
427		else
428	<	!! See if we need to update neighbor lists
428	>	!! See if we need to update neighbor lists for non pre-pair
429		call checkNeighborList(nlocal, q, listSkin, update_nlist)
430		endif
431	<
431	>
432	>
433	>
434	>
435	>
436	>	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
437	>
438	>
439	>
440	>
441	>
442		#ifdef IS_MPI
443
444		if (update_nlist) then
267	–
445		!! save current configuration, construct neighbor list,
446		!! and calculate forces
447		call saveNeighborList(nlocal, q)
#	Line 330 \| Line 507 \| contains
507		#else
508
509		if (update_nlist) then
510	<
510	>
511		! save current configuration, contruct neighbor list,
512		! and calculate forces
513		call saveNeighborList(natoms, q)
#	Line 398 \| Line 575 \| contains
575		#endif
576
577		! phew, done with main loop.
578	<
578	>
579	>	!! Do timing
580	>	#ifdef PROFILE
581	>	call cpu_time(forceTimeFinal)
582	>	forceTime = forceTime + forceTimeFinal - forceTimeInitial
583	>	#endif
584	>
585	>
586		#ifdef IS_MPI
587		!!distribute forces
588
#	Line 512 \| Line 696 \| contains
696		endif
697
698		#endif
699	<
699	>
700	>
701	>	endif
702	>
703	>	#endif
704	>
705		end subroutine do_force_loop
706
707		subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
#	Line 531 \| Line 720 \| contains
720		logical :: is_LJ_i, is_LJ_j
721		logical :: is_DP_i, is_DP_j
722		logical :: is_GB_i, is_GB_j
723	+	logical :: is_EAM_i,is_EAM_j
724		logical :: is_Sticky_i, is_Sticky_j
725		integer :: me_i, me_j
726
#	Line 564 \| Line 754 \| contains
754		call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
755
756		if ( is_DP_i .and. is_DP_j ) then
567	–
757		call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
758		do_pot, do_stress)
759		if (FF_uses_RF .and. SimUsesRF()) then
#	Line 589 \| Line 778 \| contains
778
779		if (FF_uses_GB .and. SimUsesGB()) then
780
781	+
782		call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
783		call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
784
#	Line 599 \| Line 789 \| contains
789		endif
790
791
792	+
793	+	if (FF_uses_EAM .and. SimUsesEAM()) then
794	+	call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
795	+	call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
796	+
797	+	if ( is_EAM_i .and. is_EAM_j ) &
798	+	call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
799	+	endif
800	+
801
802	+
803	+
804		end subroutine do_pair
805
806
807	+
808	+	subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
809	+	real( kind = dp ) :: pot
810	+	real( kind = dp ), dimension(3,getNlocal()) :: u_l
811	+	real (kind=dp), dimension(9,getNlocal()) :: A
812	+	real (kind=dp), dimension(3,getNlocal()) :: f
813	+	real (kind=dp), dimension(3,getNlocal()) :: t
814	+
815	+	logical, intent(inout) :: do_pot, do_stress
816	+	integer, intent(in) :: i, j
817	+	real ( kind = dp ), intent(inout) :: rijsq
818	+	real ( kind = dp ) :: r
819	+	real ( kind = dp ), intent(inout) :: d(3)
820	+
821	+	logical :: is_EAM_i, is_EAM_j
822	+
823	+	integer :: me_i, me_j
824	+
825	+	r = sqrt(rijsq)
826	+
827	+
828	+	#ifdef IS_MPI
829	+	if (tagRow(i) .eq. tagColumn(j)) then
830	+	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
831	+	endif
832	+
833	+	me_i = atid_row(i)
834	+	me_j = atid_col(j)
835	+
836	+	#else
837	+
838	+	me_i = atid(i)
839	+	me_j = atid(j)
840	+
841	+	#endif
842	+
843	+	if (FF_uses_EAM .and. SimUsesEAM()) then
844	+	call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
845	+	call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
846	+
847	+	if ( is_EAM_i .and. is_EAM_j ) &
848	+	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
849	+	endif
850	+
851	+	end subroutine do_prepair
852	+
853	+
854	+
855	+
856	+	subroutine do_preforce(nlocal,pot)
857	+	integer :: nlocal
858	+	real( kind = dp ) :: pot
859	+
860	+	if (FF_uses_EAM .and. SimUsesEAM()) then
861	+	call calc_EAM_preforce_Frho(nlocal,pot)
862	+	endif
863	+
864	+
865	+	end subroutine do_preforce
866	+
867	+
868		subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
869
870		real (kind = dp), dimension(3) :: q_i
#	Line 660 \| Line 922 \| contains
922		error = 0
923		! Make sure we are properly initialized.
924		if (.not. do_forces_initialized) then
925	+	write(,) "Forces not initialized"
926		error = -1
927		return
928		endif
#	Line 707 \| Line 970 \| contains
970
971		#endif
972
973	+
974	+	if (FF_uses_EAM .and. SimUsesEAM()) then
975	+	call clean_EAM()
976	+	endif
977	+
978	+
979	+
980	+
981	+
982		rf = 0.0_dp
983		tau_Temp = 0.0_dp
984		virial_Temp = 0.0_dp
#	Line 819 \| Line 1091 \| end module do_Forces
1091		doesit = FF_uses_RF
1092		end function FF_RequiresPostpairCalc
1093
1094	+	#ifdef PROFILE
1095	+	function getforcetime() return(totalforcetime)
1096	+	real(kind=dp) :: totalforcetime
1097	+	totalforcetime = forcetime
1098	+	end function getforcetime
1099	+	#endif
1100	+
1101	+	!! This cleans componets of force arrays belonging only to fortran
1102	+
1103		end module do_Forces

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents): Revision 626 by mmeineke, Wed Jul 16 21:30:56 2003 UTC vs. Revision 883 by chuckv, Thu Dec 18 20:46:45 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 626 by mmeineke, Wed Jul 16 21:30:56 2003 UTC vs.
Revision 883 by chuckv, Thu Dec 18 20:46:45 2003 UTC