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

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 470 by chuckv, Mon Apr 7 20:50:46 2003 UTC vs.
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.11 2003-04-07 20:50:46 chuckv Exp $, $Date: 2003-04-07 20:50:46 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $
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
#	Line 171 \| Line 218 \| contains
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
#	Line 197 \| Line 247 \| contains
247		nlocal = getNlocal()
248		natoms = nlocal
249		#endif
250	<
201	<	call getRcut(rcut,rc2=rcutsq)
202	<	call getRlist(rlist,rlistsq)
203	<
250	>
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 258 \| contains
258
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 228 \| 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(nlocal, q)
#	Line 253 \| Line 525 \| contains
525		nlist = 0
526
527		do i = 1, nrow
528	+
529		point(i) = nlist + 1
530
531		inner: do j = 1, ncol
#	Line 261 \| 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 277 \| 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_local)
283	<	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 311 \| 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(natoms, q)
#	Line 330 \| Line 602 \| contains
602		call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
603
604
605	<	if (rijsq < rlistsq) then
605	>	if (rijsq < rlistsq) then
606
607		nlist = nlist + 1
608
#	Line 346 \| Line 618 \| contains
618
619		list(nlist) = j
620
621	<	if (rijsq < rcutsq) then
350	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
621	>	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
622		u_l, A, f, t, pot)
623	<	endif
623	>
624		endif
625		enddo inner
626		enddo
#	Line 380 \| 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
#	Line 480 \| Line 758 \| contains
758		endif
759
760		if (do_stress) then
761	<	call mpi_allreduce(tau_Temp, tau,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_Temp, virial,1,mpi_double_precision,mpi_sum, &
764		mpi_comm_world,mpi_err)
#	Line 492 \| 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)
#	Line 513 \| 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 529 \| Line 814 \| contains
814		me_j = atid(j)
815
816		#endif
817	<
817	>
818		if (FF_uses_LJ .and. SimUsesLJ()) then
534	–	call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
535	–	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
542	–	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
543	–	call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
826
827	<	if ( is_DP_i .and. is_DP_j ) then
546	<
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)
830		if (FF_uses_RF .and. SimUsesRF()) then
#	Line 556 \| Line 837 \| contains
837
838		if (FF_uses_Sticky .and. SimUsesSticky()) then
839
840	<	call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
560	<	call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
561	<
562	<	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 567 \| Line 845 \| contains
845
846
847		if (FF_uses_GB .and. SimUsesGB()) then
570	–
571	–	call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
572	–	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	<	d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
944	<	int(abs(d(1:3)/box(1:3)) + 0.5_dp)
943	>	if( .not.boxIsOrthorhombic ) then
944	>	! calc the scaled coordinates.
945	>
946	>	scaled = matmul(HmatInv, d)
947	>
948	>	! wrap the scaled coordinates
949	>
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 653 \| 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
#	Line 765 \| 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 trunk/OOPSE/libmdtools/do_Forces.F90 (file contents): Revision 470 by chuckv, Mon Apr 7 20:50:46 2003 UTC vs. Revision 895 by chuckv, Mon Jan 5 22:12:11 2004 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 470 by chuckv, Mon Apr 7 20:50:46 2003 UTC vs.
Revision 895 by chuckv, Mon Jan 5 22:12:11 2004 UTC