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

Comparing:
branches/mmeineke/OOPSE/libmdtools/do_Forces.F90 (file contents), Revision 377 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
trunk/OOPSE/libmdtools/do_Forces.F90 (file contents), Revision 895 by chuckv, Mon Jan 5 22:12:11 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.1.1.1 2003-03-21 17:42:12 mmeineke Exp $, $Date: 2003-03-21 17:42:12 $, $Name: not supported by cvs2svn $, $Revision: 1.1.1.1 $
7	>	!! @version $Id: do_Forces.F90,v 1.42 2004-01-05 22:12:11 chuckv Exp $, $Date: 2004-01-05 22:12:11 $, $Name: not supported by cvs2svn $, $Revision: 1.42 $
8
9		module do_Forces
10		use force_globals
#	Line 17 \| Line 17 \| module do_Forces
17		use dipole_dipole
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 27 \| Line 30 \| module do_Forces
30		#define __FORTRAN90
31		#include "fForceField.h"
32
33	<	logical, save :: do_forces_initialized = .false.
33	>	logical, save :: do_forces_initialized = .false., haveRlist = .false.
34	>	logical, save :: havePolicies = .false.
35		logical, save :: FF_uses_LJ
36		logical, save :: FF_uses_sticky
37		logical, save :: FF_uses_dipoles
#	Line 35 \| Line 39 \| module do_Forces
39		logical, save :: FF_uses_GB
40		logical, save :: FF_uses_EAM
41
42	+	real(kind=dp), save :: rlist, rlistsq
43	+
44		public :: init_FF
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	+	integer :: nLoops
53	+	#endif
54	+
55	+	logical, allocatable :: propertyMapI(:,:)
56	+	logical, allocatable :: propertyMapJ(:,:)
57	+
58		contains
59
60	+	subroutine setRlistDF( this_rlist )
61	+
62	+	real(kind=dp) :: this_rlist
63	+
64	+	rlist = this_rlist
65	+	rlistsq = rlist * rlist
66	+
67	+	haveRlist = .true.
68	+	if( havePolicies ) do_forces_initialized = .true.
69	+
70	+	end subroutine setRlistDF
71	+
72		subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
73
74		integer, intent(in) :: LJMIXPOLICY
#	Line 87 \| Line 116 \| contains
116		!! check to make sure the FF_uses_RF setting makes sense
117
118		if (FF_uses_dipoles) then
90	–	rrf = getRrf()
91	–	rt = getRt()
92	–	call initialize_dipole(rrf, rt)
119		if (FF_uses_RF) then
120		dielect = getDielect()
121	<	call initialize_rf(rrf, rt, dielect)
121	>	call initialize_rf(dielect)
122		endif
123		else
124		if (FF_uses_RF) then
#	Line 100 \| Line 126 \| contains
126		thisStat = -1
127		return
128		endif
129	<	endif
129	>	endif
130
131		if (FF_uses_LJ) then
132
107	–	call getRcut(rcut)
108	–
133		select case (LJMIXPOLICY)
134		case (LB_MIXING_RULE)
135	<	call init_lj_FF(LB_MIXING_RULE, rcut, my_status)
135	>	call init_lj_FF(LB_MIXING_RULE, my_status)
136		case (EXPLICIT_MIXING_RULE)
137	<	call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
137	>	call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
138		case default
139		write(default_error,*) 'unknown LJ Mixing Policy!'
140		thisStat = -1
#	Line 129 \| Line 153 \| contains
153		return
154		end if
155		endif
156	+
157	+
158	+	if (FF_uses_EAM) then
159	+	call init_EAM_FF(my_status)
160	+	if (my_status /= 0) then
161	+	write(,) "init_EAM_FF returned a bad status"
162	+	thisStat = -1
163	+	return
164	+	end if
165	+	endif
166	+
167	+
168
169		if (FF_uses_GB) then
170		call check_gb_pair_FF(my_status)
#	Line 140 \| Line 176 \| contains
176
177		if (FF_uses_GB .and. FF_uses_LJ) then
178		endif
179	+	if (.not. do_forces_initialized) then
180	+	!! Create neighbor lists
181	+	call expandNeighborList(getNlocal(), my_status)
182	+	if (my_Status /= 0) then
183	+	write(default_error,*) "SimSetup: ExpandNeighborList returned error."
184	+	thisStat = -1
185	+	return
186	+	endif
187	+	endif
188	+
189
190	+	havePolicies = .true.
191	+	if( haveRlist ) do_forces_initialized = .true.
192
145	–	do_forces_initialized = .true.
146	–
193		end subroutine init_FF
194
195
#	Line 161 \| Line 207 \| contains
207		real ( kind = dp ), dimension(3,getNlocal()) :: f
208		!! Torsion array provided by C, dimensioned by getNlocal
209		real( kind = dp ), dimension(3,getNlocal()) :: t
210	+
211		!! Stress Tensor
212		real( kind = dp), dimension(9) :: tau
213		real ( kind = dp ) :: pot
214		logical ( kind = 2) :: do_pot_c, do_stress_c
215		logical :: do_pot
216		logical :: do_stress
217	<	#ifdef IS_MPI
217	>	#ifdef IS_MPI
218		real( kind = DP ) :: pot_local
219		integer :: nrow
220		integer :: ncol
221	+	integer :: nprocs
222		#endif
223		integer :: nlocal
224		integer :: natoms
225		logical :: update_nlist
226		integer :: i, j, jbeg, jend, jnab
227		integer :: nlist
228	<	real( kind = DP ) :: rijsq, rlistsq, rcutsq, rlist, rcut
228	>	real( kind = DP ) :: rijsq
229		real(kind=dp),dimension(3) :: d
230		real(kind=dp) :: rfpot, mu_i, virial
231		integer :: me_i
#	Line 186 \| Line 234 \| contains
234		integer :: listerror, error
235		integer :: localError
236
237	+	real(kind=dp) :: listSkin = 1.0
238	+
239		!! initialize local variables
240
241		#ifdef IS_MPI
242	+	pot_local = 0.0_dp
243		nlocal = getNlocal()
244		nrow = getNrow(plan_row)
245		ncol = getNcol(plan_col)
#	Line 197 \| Line 248 \| contains
248		natoms = nlocal
249		#endif
250
200	–	call getRcut(rcut,rc2=rcutsq)
201	–	call getRlist(rlist,rlistsq)
202	–
251		call check_initialization(localError)
252		if ( localError .ne. 0 ) then
253	+	call handleError("do_force_loop","Not Initialized")
254		error = -1
255		return
256		end if
#	Line 210 \| Line 259 \| contains
259		do_pot = do_pot_c
260		do_stress = do_stress_c
261
262	+
263	+	#ifdef IS_MPI
264	+	if (.not.allocated(propertyMapI)) then
265	+	allocate(propertyMapI(5,getNrow())
266	+	endif
267	+
268	+	do i = 1, nrow
269	+	me_i = atid_row(i)
270	+	#else
271	+	if (.not.allocated(propertyMapI)) then
272	+	allocate(propertyMapI(5,getNlocal())
273	+	endif
274	+
275	+	do i = 1, natoms
276	+	me_i = atid(i)
277	+	#endif
278	+
279	+	propertyMapI(1:5,i) = .false.
280	+
281	+	call getElementProperty(atypes, me_i, "propertyPack", propPack_i)
282	+
283	+	! unpack the properties
284	+
285	+	if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
286	+	propertyMapI(1, i) = .true.
287	+	if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
288	+	propertyMapI(2, i) = .true.
289	+	if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
290	+	propertyMapI(3, i) = .true.
291	+	if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
292	+	propertyMapI(4, i) = .true.
293	+	if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
294	+	propertyMapI(5, i) = .true.
295	+
296	+	end do
297	+
298	+	#ifdef IS_MPI
299	+	if (.not.allocated(propertyMapJ)) then
300	+	allocate(propertyMapJ(5,getNcol())
301	+	endif
302	+
303	+	do j = 1, ncol
304	+	me_j = atid_col(j)
305	+	#else
306	+	if (.not.allocated(propertyMapJ)) then
307	+	allocate(propertyMapJ(5,getNlocal())
308	+	endif
309	+
310	+	do j = 1, natoms
311	+	me_j = atid(j)
312	+	#endif
313	+
314	+	propertyMapJ(1:5,j) = .false.
315	+
316	+	call getElementProperty(atypes, me_j, "propertyPack", propPack_j)
317	+
318	+	! unpack the properties
319	+
320	+	if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
321	+	propertyMapJ(1, j) = .true.
322	+	if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
323	+	propertyMapJ(2, j) = .true.
324	+	if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
325	+	propertyMapJ(3, j) = .true.
326	+	if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
327	+	propertyMapJ(4, j) = .true.
328	+	if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
329	+	propertyMapJ(5, j) = .true.
330	+
331	+	end do
332	+
333		! Gather all information needed by all force loops:
334
335		#ifdef IS_MPI
#	Line 226 \| Line 346 \| contains
346		endif
347
348		#endif
349	<
349	>
350	>	!! Begin force loop timing:
351	>	#ifdef PROFILE
352	>	call cpu_time(forceTimeInitial)
353	>	nloops = nloops + 1
354	>	#endif
355	>
356		if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
357		!! See if we need to update neighbor lists
358	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
358	>	call checkNeighborList(nlocal, q, listSkin, update_nlist)
359		!! if_mpi_gather_stuff_for_prepair
360		!! do_prepair_loop_if_needed
361		!! if_mpi_scatter_stuff_from_prepair
362		!! if_mpi_gather_stuff_from_prepair_to_main_loop
363	<	else
364	<	!! See if we need to update neighbor lists
365	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
363	>
364	>	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
365	>	#ifdef IS_MPI
366	>
367	>	if (update_nlist) then
368	>
369	>	!! save current configuration, construct neighbor list,
370	>	!! and calculate forces
371	>	call saveNeighborList(nlocal, q)
372	>
373	>	neighborListSize = size(list)
374	>	nlist = 0
375	>
376	>	do i = 1, nrow
377	>	point(i) = nlist + 1
378	>
379	>	prepair_inner: do j = 1, ncol
380	>
381	>	if (skipThisPair(i,j)) cycle prepair_inner
382	>
383	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
384	>
385	>	if (rijsq < rlistsq) then
386	>
387	>	nlist = nlist + 1
388	>
389	>	if (nlist > neighborListSize) then
390	>	call expandNeighborList(nlocal, listerror)
391	>	if (listerror /= 0) then
392	>	error = -1
393	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
394	>	return
395	>	end if
396	>	neighborListSize = size(list)
397	>	endif
398	>
399	>	list(nlist) = j
400	>	call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)
401	>	endif
402	>	enddo prepair_inner
403	>	enddo
404	>
405	>	point(nrow + 1) = nlist + 1
406	>
407	>	else !! (of update_check)
408	>
409	>	! use the list to find the neighbors
410	>	do i = 1, nrow
411	>	JBEG = POINT(i)
412	>	JEND = POINT(i+1) - 1
413	>	! check thiat molecule i has neighbors
414	>	if (jbeg .le. jend) then
415	>
416	>	do jnab = jbeg, jend
417	>	j = list(jnab)
418	>
419	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
420	>	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
421	>	u_l, A, f, t, pot_local)
422	>
423	>	enddo
424	>	endif
425	>	enddo
426		endif
427
428	<	#ifdef IS_MPI
428	>	#else
429
430		if (update_nlist) then
431
432	+	! save current configuration, contruct neighbor list,
433	+	! and calculate forces
434	+	call saveNeighborList(natoms, q)
435	+
436	+	neighborListSize = size(list)
437	+
438	+	nlist = 0
439	+
440	+	do i = 1, natoms-1
441	+	point(i) = nlist + 1
442	+
443	+	prepair_inner: do j = i+1, natoms
444	+
445	+	if (skipThisPair(i,j)) cycle prepair_inner
446	+
447	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
448	+
449	+
450	+	if (rijsq < rlistsq) then
451	+
452	+
453	+	nlist = nlist + 1
454	+
455	+	if (nlist > neighborListSize) then
456	+	call expandNeighborList(natoms, listerror)
457	+	if (listerror /= 0) then
458	+	error = -1
459	+	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
460	+	return
461	+	end if
462	+	neighborListSize = size(list)
463	+	endif
464	+
465	+	list(nlist) = j
466	+
467	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
468	+	u_l, A, f, t, pot)
469	+
470	+	endif
471	+	enddo prepair_inner
472	+	enddo
473	+
474	+	point(natoms) = nlist + 1
475	+
476	+	else !! (update)
477	+
478	+	! use the list to find the neighbors
479	+	do i = 1, natoms-1
480	+	JBEG = POINT(i)
481	+	JEND = POINT(i+1) - 1
482	+	! check thiat molecule i has neighbors
483	+	if (jbeg .le. jend) then
484	+
485	+	do jnab = jbeg, jend
486	+	j = list(jnab)
487	+
488	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
489	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
490	+	u_l, A, f, t, pot)
491	+
492	+	enddo
493	+	endif
494	+	enddo
495	+	endif
496	+	#endif
497	+	!! Do rest of preforce calculations
498	+	!! do necessary preforce calculations
499	+	call do_preforce(nlocal,pot)
500	+	! we have already updated the neighbor list set it to false...
501	+	update_nlist = .false.
502	+	else
503	+	!! See if we need to update neighbor lists for non pre-pair
504	+	call checkNeighborList(nlocal, q, listSkin, update_nlist)
505	+	endif
506	+
507	+
508	+
509	+
510	+
511	+	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
512	+
513	+
514	+
515	+
516	+
517	+	#ifdef IS_MPI
518	+
519	+	if (update_nlist) then
520		!! save current configuration, construct neighbor list,
521		!! and calculate forces
522	<	call saveNeighborList(q)
522	>	call saveNeighborList(nlocal, q)
523
524		neighborListSize = size(list)
525		nlist = 0
526
527		do i = 1, nrow
528	+
529		point(i) = nlist + 1
530
531		inner: do j = 1, ncol
#	Line 259 \| Line 534 \| contains
534
535		call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
536
537	<	if (rijsq < rlistsq) then
537	>	if (rijsq < rlistsq) then
538
539		nlist = nlist + 1
540
#	Line 275 \| Line 550 \| contains
550
551		list(nlist) = j
552
553	<	if (rijsq < rcutsq) then
554	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
555	<	u_l, A, f, t,pot)
281	<	endif
553	>	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
554	>	u_l, A, f, t, pot_local)
555	>
556		endif
557		enddo inner
558		enddo
#	Line 299 \| Line 573 \| contains
573
574		call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
575		call do_pair(i, j, rijsq, d, do_pot, do_stress, &
576	<	u_l, A, f, t,pot)
576	>	u_l, A, f, t, pot_local)
577
578		enddo
579		endif
#	Line 309 \| Line 583 \| contains
583		#else
584
585		if (update_nlist) then
586	<
586	>
587		! save current configuration, contruct neighbor list,
588		! and calculate forces
589	<	call saveNeighborList(q)
589	>	call saveNeighborList(natoms, q)
590
591		neighborListSize = size(list)
592
#	Line 323 \| Line 597 \| contains
597
598		inner: do j = i+1, natoms
599
600	<	if (skipThisPair(i,j)) cycle inner
601	<
600	>	if (skipThisPair(i,j)) cycle inner
601	>
602		call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
603
604	<	if (rijsq < rlistsq) then
604	>
605	>	if (rijsq < rlistsq) then
606
607		nlist = nlist + 1
608
#	Line 343 \| Line 618 \| contains
618
619		list(nlist) = j
620
621	<	if (rijsq < rcutsq) then
622	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
623	<	u_l, A, f, t,pot)
349	<	endif
621	>	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
622	>	u_l, A, f, t, pot)
623	>
624		endif
625		enddo inner
626		enddo
#	Line 367 \| Line 641 \| contains
641
642		call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
643		call do_pair(i, j, rijsq, d, do_pot, do_stress, &
644	<	u_l, A, f, t,pot)
644	>	u_l, A, f, t, pot)
645
646		enddo
647		endif
#	Line 377 \| Line 651 \| contains
651		#endif
652
653		! phew, done with main loop.
654	<
654	>
655	>	!! Do timing
656	>	#ifdef PROFILE
657	>	call cpu_time(forceTimeFinal)
658	>	forceTime = forceTime + forceTimeFinal - forceTimeInitial
659	>	#endif
660	>
661	>
662		#ifdef IS_MPI
663		!!distribute forces
664	<
665	<	call scatter(f_Row,f,plan_row3d)
664	>
665	>	f_temp = 0.0_dp
666	>	call scatter(f_Row,f_temp,plan_row3d)
667	>	do i = 1,nlocal
668	>	f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
669	>	end do
670	>
671	>	f_temp = 0.0_dp
672		call scatter(f_Col,f_temp,plan_col3d)
673		do i = 1,nlocal
674		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
675		end do
676
677		if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
678	<	call scatter(t_Row,t,plan_row3d)
678	>	t_temp = 0.0_dp
679	>	call scatter(t_Row,t_temp,plan_row3d)
680	>	do i = 1,nlocal
681	>	t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
682	>	end do
683	>	t_temp = 0.0_dp
684		call scatter(t_Col,t_temp,plan_col3d)
685
686		do i = 1,nlocal
#	Line 399 \| Line 691 \| contains
691		if (do_pot) then
692		! scatter/gather pot_row into the members of my column
693		call scatter(pot_Row, pot_Temp, plan_row)
694	<
694	>
695		! scatter/gather pot_local into all other procs
696		! add resultant to get total pot
697		do i = 1, nlocal
698		pot_local = pot_local + pot_Temp(i)
699		enddo
700	+
701	+	pot_Temp = 0.0_DP
702
409	–	pot_Temp = 0.0_DP
410	–
703		call scatter(pot_Col, pot_Temp, plan_col)
704		do i = 1, nlocal
705		pot_local = pot_local + pot_Temp(i)
706		enddo
707	<
707	>
708		endif
709		#endif
710
#	Line 460 \| Line 752 \| contains
752		#ifdef IS_MPI
753
754		if (do_pot) then
755	<	pot = pot_local
755	>	pot = pot + pot_local
756		!! we assume the c code will do the allreduce to get the total potential
757		!! we could do it right here if we needed to...
758		endif
759
760		if (do_stress) then
761	<	call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
761	>	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
762		mpi_comm_world,mpi_err)
763	<	call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
763	>	call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
764		mpi_comm_world,mpi_err)
765		endif
766
#	Line 478 \| Line 770 \| contains
770		tau = tau_Temp
771		virial = virial_Temp
772		endif
773	<
773	>
774		#endif
775
776	+
777	+
778		end subroutine do_force_loop
779
780	<	subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t,pot)
780	>	subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
781
782		real( kind = dp ) :: pot
783	<	real( kind = dp ), dimension(:,:) :: u_l
784	<	real (kind=dp), dimension(:,:) :: A
785	<	real (kind=dp), dimension(:,:) :: f
786	<	real (kind=dp), dimension(:,:) :: t
783	>	real( kind = dp ), dimension(3,getNlocal()) :: u_l
784	>	real (kind=dp), dimension(9,getNlocal()) :: A
785	>	real (kind=dp), dimension(3,getNlocal()) :: f
786	>	real (kind=dp), dimension(3,getNlocal()) :: t
787
788		logical, intent(inout) :: do_pot, do_stress
789		integer, intent(in) :: i, j
#	Line 499 \| Line 793 \| contains
793		logical :: is_LJ_i, is_LJ_j
794		logical :: is_DP_i, is_DP_j
795		logical :: is_GB_i, is_GB_j
796	+	logical :: is_EAM_i,is_EAM_j
797		logical :: is_Sticky_i, is_Sticky_j
798		integer :: me_i, me_j
799	<
799	>	integer :: propPack_i
800	>	integer :: propPack_j
801		r = sqrt(rijsq)
802
803		#ifdef IS_MPI
804	+	if (tagRow(i) .eq. tagColumn(j)) then
805	+	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
806	+	endif
807
808		me_i = atid_row(i)
809		me_j = atid_col(j)
#	Line 515 \| Line 814 \| contains
814		me_j = atid(j)
815
816		#endif
817	<
817	>
818		if (FF_uses_LJ .and. SimUsesLJ()) then
520	–	call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
521	–	call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
819
820	<	if ( is_LJ_i .and. is_LJ_j ) &
820	>	if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) &
821		call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
822	+
823		endif
824
825		if (FF_uses_dipoles .and. SimUsesDipoles()) then
528	–	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
529	–	call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
826
827	<	if ( is_DP_i .and. is_DP_j ) then
532	<
827	>	if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then
828		call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
829		do_pot, do_stress)
535	–
830		if (FF_uses_RF .and. SimUsesRF()) then
537	–
831		call accumulate_rf(i, j, r, u_l)
832		call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
540	–
833		endif
834
835		endif
#	Line 545 \| Line 837 \| contains
837
838		if (FF_uses_Sticky .and. SimUsesSticky()) then
839
840	<	call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
549	<	call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
550	<
551	<	if ( is_Sticky_i .and. is_Sticky_j ) then
840	>	if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then
841		call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
842		do_pot, do_stress)
843		endif
#	Line 556 \| Line 845 \| contains
845
846
847		if (FF_uses_GB .and. SimUsesGB()) then
559	–
560	–	call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
561	–	call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
848
849	<	if ( is_GB_i .and. is_GB_j ) then
849	>	if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then
850		call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
851		do_pot, do_stress)
852		endif
853	+
854		endif
855
856	+
857	+
858	+	if (FF_uses_EAM .and. SimUsesEAM()) then
859	+
860	+	if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then
861	+	call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
862	+	endif
863	+
864	+	endif
865	+
866		end subroutine do_pair
867
868
869	+
870	+	subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
871	+	real( kind = dp ) :: pot
872	+	real( kind = dp ), dimension(3,getNlocal()) :: u_l
873	+	real (kind=dp), dimension(9,getNlocal()) :: A
874	+	real (kind=dp), dimension(3,getNlocal()) :: f
875	+	real (kind=dp), dimension(3,getNlocal()) :: t
876	+
877	+	logical, intent(inout) :: do_pot, do_stress
878	+	integer, intent(in) :: i, j
879	+	real ( kind = dp ), intent(inout) :: rijsq
880	+	real ( kind = dp ) :: r
881	+	real ( kind = dp ), intent(inout) :: d(3)
882	+
883	+	logical :: is_EAM_i, is_EAM_j
884	+
885	+	integer :: me_i, me_j
886	+
887	+	r = sqrt(rijsq)
888	+
889	+
890	+	#ifdef IS_MPI
891	+	if (tagRow(i) .eq. tagColumn(j)) then
892	+	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
893	+	endif
894	+
895	+	me_i = atid_row(i)
896	+	me_j = atid_col(j)
897	+
898	+	#else
899	+
900	+	me_i = atid(i)
901	+	me_j = atid(j)
902	+
903	+	#endif
904	+
905	+	if (FF_uses_EAM .and. SimUsesEAM()) then
906	+	call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
907	+	call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
908	+
909	+	if ( is_EAM_i .and. is_EAM_j ) &
910	+	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
911	+	endif
912	+
913	+	end subroutine do_prepair
914	+
915	+
916	+
917	+
918	+	subroutine do_preforce(nlocal,pot)
919	+	integer :: nlocal
920	+	real( kind = dp ) :: pot
921	+
922	+	if (FF_uses_EAM .and. SimUsesEAM()) then
923	+	call calc_EAM_preforce_Frho(nlocal,pot)
924	+	endif
925	+
926	+
927	+	end subroutine do_preforce
928	+
929	+
930		subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
931
932		real (kind = dp), dimension(3) :: q_i
933		real (kind = dp), dimension(3) :: q_j
934		real ( kind = dp ), intent(out) :: r_sq
935	<	real( kind = dp ) :: d(3)
936	<	real( kind = dp ) :: d_old(3)
937	<	d(1:3) = q_i(1:3) - q_j(1:3)
938	<	d_old = d
935	>	real( kind = dp ) :: d(3), scaled(3)
936	>	integer i
937	>
938	>	d(1:3) = q_j(1:3) - q_i(1:3)
939	>
940		! Wrap back into periodic box if necessary
941		if ( SimUsesPBC() ) then
942	+
943	+	if( .not.boxIsOrthorhombic ) then
944	+	! calc the scaled coordinates.
945	+
946	+	scaled = matmul(HmatInv, d)
947	+
948	+	! wrap the scaled coordinates
949
950	<	d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
951	<	int(abs(d(1:3)/box(1:3)) + 0.5_dp)
950	>	scaled = scaled - anint(scaled)
951	>
952
953	+	! calc the wrapped real coordinates from the wrapped scaled
954	+	! coordinates
955	+
956	+	d = matmul(Hmat,scaled)
957	+
958	+	else
959	+	! calc the scaled coordinates.
960	+
961	+	do i = 1, 3
962	+	scaled(i) = d(i) * HmatInv(i,i)
963	+
964	+	! wrap the scaled coordinates
965	+
966	+	scaled(i) = scaled(i) - anint(scaled(i))
967	+
968	+	! calc the wrapped real coordinates from the wrapped scaled
969	+	! coordinates
970	+
971	+	d(i) = scaled(i)*Hmat(i,i)
972	+	enddo
973	+	endif
974	+
975		endif
976	+
977		r_sq = dot_product(d,d)
978	<
978	>
979		end subroutine get_interatomic_vector
980	<
980	>
981		subroutine check_initialization(error)
982		integer, intent(out) :: error
983
984		error = 0
985		! Make sure we are properly initialized.
986		if (.not. do_forces_initialized) then
987	+	write(,) "Forces not initialized"
988		error = -1
989		return
990		endif
#	Line 642 \| Line 1032 \| contains
1032
1033		#endif
1034
1035	+
1036	+	if (FF_uses_EAM .and. SimUsesEAM()) then
1037	+	call clean_EAM()
1038	+	endif
1039	+
1040	+
1041	+
1042	+
1043	+
1044		rf = 0.0_dp
1045		tau_Temp = 0.0_dp
1046		virial_Temp = 0.0_dp
648	–
1047		end subroutine zero_work_arrays
1048
1049		function skipThisPair(atom1, atom2) result(skip_it)
652	–
1050		integer, intent(in) :: atom1
1051		integer, intent(in), optional :: atom2
1052		logical :: skip_it
1053		integer :: unique_id_1, unique_id_2
1054	+	integer :: me_i,me_j
1055		integer :: i
1056
1057		skip_it = .false.
#	Line 671 \| Line 1069 \| contains
1069		!! in the normal loop, the atom numbers are unique
1070		unique_id_1 = atom1
1071		#endif
1072	<
1072	>
1073		!! We were called with only one atom, so just check the global exclude
1074		!! list for this atom
1075		if (.not. present(atom2)) then
#	Line 689 \| Line 1087 \| contains
1087		#else
1088		unique_id_2 = atom2
1089		#endif
1090	<
1090	>
1091		#ifdef IS_MPI
1092		!! this situation should only arise in MPI simulations
1093		if (unique_id_1 == unique_id_2) then
#	Line 699 \| Line 1097 \| contains
1097
1098		!! this prevents us from doing the pair on multiple processors
1099		if (unique_id_1 < unique_id_2) then
1100	<	if (mod(unique_id_1 + unique_id_2,2) == 0) skip_it = .true.
1101	<	return
1100	>	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1101	>	skip_it = .true.
1102	>	return
1103	>	endif
1104		else
1105	<	if (mod(unique_id_1 + unique_id_2,2) == 1) skip_it = .true.
1106	<	return
1105	>	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1106	>	skip_it = .true.
1107	>	return
1108	>	endif
1109		endif
1110		#endif
1111	<
1111	>
1112		!! the rest of these situations can happen in all simulations:
1113		do i = 1, nExcludes_global
1114		if ((excludesGlobal(i) == unique_id_1) .or. &
#	Line 715 \| Line 1117 \| contains
1117		return
1118		endif
1119		enddo
1120	<
1120	>
1121		do i = 1, nExcludes_local
1122		if (excludesLocal(1,i) == unique_id_1) then
1123		if (excludesLocal(2,i) == unique_id_2) then
#	Line 751 \| Line 1153 \| end module do_Forces
1153		doesit = FF_uses_RF
1154		end function FF_RequiresPostpairCalc
1155
1156	+	#ifdef PROFILE
1157	+	function getforcetime() result(totalforcetime)
1158	+	real(kind=dp) :: totalforcetime
1159	+	totalforcetime = forcetime
1160	+	end function getforcetime
1161	+	#endif
1162	+
1163	+	!! This cleans componets of force arrays belonging only to fortran
1164	+
1165		end module do_Forces

Diff Legend

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

Comparing: branches/mmeineke/OOPSE/libmdtools/do_Forces.F90 (file contents), Revision 377 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs. trunk/OOPSE/libmdtools/do_Forces.F90 (file contents), Revision 895 by chuckv, Mon Jan 5 22:12:11 2004 UTC

Diff Legend

Comparing:
branches/mmeineke/OOPSE/libmdtools/do_Forces.F90 (file contents), Revision 377 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
trunk/OOPSE/libmdtools/do_Forces.F90 (file contents), Revision 895 by chuckv, Mon Jan 5 22:12:11 2004 UTC