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

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

Diff Legend

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

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

Diff Legend

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