[ViewVC] Diff of: group/trunk/mdtools/md_code/lj

Comparing trunk/mdtools/md_code/lj_FF.F90 (file contents):
Revision 252 by chuckv, Tue Jan 28 22:16:55 2003 UTC vs.
Revision 281 by chuckv, Mon Feb 24 21:25:15 2003 UTC

#	Line 2 \| Line 2
2		!! Corresponds to the force field defined in lj_FF.cpp
3		!! @author Charles F. Vardeman II
4		!! @author Matthew Meineke
5	<	!! @version $Id: lj_FF.F90,v 1.12 2003-01-28 22:16:55 chuckv Exp $, $Date: 2003-01-28 22:16:55 $, $Name: not supported by cvs2svn $, $Revision: 1.12 $
5	>	!! @version $Id: lj_FF.F90,v 1.19 2003-02-24 21:25:15 chuckv Exp $, $Date: 2003-02-24 21:25:15 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $
6
7
8
#	Line 36 \| Line 36 \| module lj_ff
36		integer :: nListAllocs = 0
37		integer, parameter :: maxListAllocs = 5
38
39	+	logical, save :: firstTime = .True.
40
41		!! Atype identity pointer lists
42		#ifdef IS_MPI
43		!! Row lj_atype pointer list
44	<	type (lj_atypePtr), dimension(:), pointer :: identPtrListRow => null()
44	>	type (lj_identPtrList), dimension(:), pointer :: identPtrListRow => null()
45		!! Column lj_atype pointer list
46	<	type (lj_atypePtr), dimension(:), pointer :: identPtrListColumn => null()
46	>	type (lj_identPtrList), dimension(:), pointer :: identPtrListColumn => null()
47		#else
48		type( lj_identPtrList ), dimension(:), pointer :: identPtrList => null()
49		#endif
#	Line 123 \| Line 124 \| contains
124
125		integer :: thisStat
126		integer :: i
127	+
128	+	integer :: myNode
129		#ifdef IS_MPI
130		integer, allocatable, dimension(:) :: identRow
131		integer, allocatable, dimension(:) :: identCol
132		integer :: nrow
133		integer :: ncol
131	–	integer :: alloc_stat
134		#endif
135		status = 0
136	+
137	+
138	+
139	+
140		!! if were're not in MPI, we just update ljatypePtrList
141		#ifndef IS_MPI
142		call create_IdentPtrlst(ident,ljListHead,identPtrList,thisStat)
#	Line 154 \| Line 160 \| contains
160
161		! if were're in MPI, we also have to worry about row and col lists
162		#else
163	+
164		! We can only set up forces if mpiSimulation has been setup.
165		if (.not. isMPISimSet()) then
166	+	write(default_error,*) "MPI is not set"
167		status = -1
168		return
169		endif
170	<	nrow = getNrow()
171	<	ncol = getNcol()
170	>	nrow = getNrow(plan_row)
171	>	ncol = getNcol(plan_col)
172	>	mynode = getMyNode()
173		!! Allocate temperary arrays to hold gather information
174		allocate(identRow(nrow),stat=alloc_stat)
175		if (alloc_stat /= 0 ) then
#	Line 173 \| Line 182 \| contains
182		status = -1
183		return
184		endif
185	+
186		!! Gather idents into row and column idents
187	+
188		call gather(ident,identRow,plan_row)
189		call gather(ident,identCol,plan_col)
190	<
190	>
191	>
192		!! Create row and col pointer lists
193	<	call create_IdentPtrlst(identRow,ljListTail,identPtrListRow,thisStat)
193	>
194	>	call create_IdentPtrlst(identRow,ljListHead,identPtrListRow,thisStat)
195		if (thisStat /= 0 ) then
196		status = -1
197		return
198		endif
199	<
200	<	call create_IdentPtrlst(identCol,ljListTail,identPtrListCol,thisStat)
199	>
200	>	call create_IdentPtrlst(identCol,ljListHead,identPtrListColumn,thisStat)
201		if (thisStat /= 0 ) then
202		status = -1
203		return
204		endif
205
206		!! free temporary ident arrays
207	<	deallocate(identCol)
208	<	deallocate(identRow)
207	>	if (allocated(identCol)) then
208	>	deallocate(identCol)
209	>	end if
210	>	if (allocated(identCol)) then
211	>	deallocate(identRow)
212	>	endif
213
197	–
214		!! Allocate neighbor lists for mpi simulations.
215		if (.not. allocated(point)) then
216		allocate(point(nrow),stat=alloc_stat)
#	Line 233 \| Line 249 \| contains
249		endif
250		isljFFinit = .true.
251
252	+
253		end subroutine init_ljFF
254
255
#	Line 331 \| Line 348 \| contains
348
349		!! FORCE routine Calculates Lennard Jones forces.
350		!------------------------------------------------------------->
351	<	subroutine do_lj_ff(q,f,potE,do_pot)
351	>	subroutine do_lj_ff(q,f,potE,tau,do_pot)
352		!! Position array provided by C, dimensioned by getNlocal
353		real ( kind = dp ), dimension(3,getNlocal()) :: q
354		!! Force array provided by C, dimensioned by getNlocal
355		real ( kind = dp ), dimension(3,getNlocal()) :: f
356	+	!! Stress Tensor
357	+	real( kind = dp), dimension(9) :: tau
358	+	real( kind = dp), dimension(9) :: tauTemp
359		real ( kind = dp ) :: potE
360		logical ( kind = 2) :: do_pot
361
362		type(lj_atype), pointer :: ljAtype_i
363		type(lj_atype), pointer :: ljAtype_j
364
365	<	#ifdef MPI
366	<	real( kind = DP ), dimension(3,getNcol()) :: efr
367	<	real( kind = DP ) :: pot_local
348	<	#else
349	<	real( kind = DP ), dimension(3,getNlocal()) :: efr
350	<	#endif
365	>
366	>
367	>
368
369	<	!! Local arrays needed for MPI
369	>
370		#ifdef IS_MPI
371	<	real(kind = dp), dimension(3:getNrow()) :: qRow
355	<	real(kind = dp), dimension(3:getNcol()) :: qCol
371	>	real( kind = DP ) :: pot_local
372
373	<	real(kind = dp), dimension(3:getNrow()) :: fRow
374	<	real(kind = dp), dimension(3:getNcol()) :: fCol
373	>	!! Local arrays needed for MPI
374	>	real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
375	>	real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
376
377	<	real(kind = dp), dimension(getNrow()) :: eRow
378	<	real(kind = dp), dimension(getNcol()) :: eCol
377	>	real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
378	>	real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
379	>	real(kind = dp), dimension(3,getNlocal()) :: fMPITemp = 0.0_dp
380
381	<	real(kind = dp), dimension(getNlocal()) :: eTemp
382	<	#endif
381	>	real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
382	>	real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
383
384	+	real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
385
386	+	#endif
387
388	+
389	+
390		real( kind = DP ) :: pe
391		logical :: update_nlist
392
#	Line 374 \| Line 396 \| contains
396		integer :: j_start
397		integer :: tag_i,tag_j
398		real( kind = DP ) :: r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
399	+	real( kind = dp ) :: fx,fy,fz
400	+	real( kind = DP ) :: drdx, drdy, drdz
401		real( kind = DP ) :: rxi, ryi, rzi, rxij, ryij, rzij, rijsq
402		real( kind = DP ) :: rlistsq, rcutsq,rlist,rcut
403
#	Line 386 \| Line 410 \| contains
410		integer :: nrow
411		integer :: ncol
412		integer :: natoms
413	+	!! should we calculate the stress tensor
414	+	logical :: do_stress = .false.
415	+
416
417	+
418		! Make sure we are properly initialized.
419		if (.not. isljFFInit) then
420		write(default_error,*) "ERROR: lj_FF has not been properly initialized"
#	Line 403 \| Line 431 \| contains
431		natoms = getNlocal()
432		call getRcut(rcut,rcut2=rcutsq)
433		call getRlist(rlist,rlistsq)
434	+	!! Find ensemble
435	+	if (isEnsemble("NPT")) do_stress = .true.
436	+
437		#ifndef IS_MPI
438		nrow = natoms - 1
439		ncol = natoms
#	Line 416 \| Line 447 \| contains
447
448		!! See if we need to update neighbor lists
449		call check(q,update_nlist)
450	+	if (firstTime) then
451	+	update_nlist = .true.
452	+	firstTime = .false.
453	+	endif
454
455		!--------------WARNING...........................
456		! Zero variables, NOTE:::: Forces are zeroed in C
#	Line 427 \| Line 462 \| contains
462		pe = 0.0E0_DP
463
464		#else
465	<	f_row = 0.0E0_DP
466	<	f_col = 0.0E0_DP
465	>	fRow = 0.0E0_DP
466	>	fCol = 0.0E0_DP
467
468		pe_local = 0.0E0_DP
469
#	Line 436 \| Line 471 \| contains
471		eCol = 0.0E0_DP
472		eTemp = 0.0E0_DP
473		#endif
439	–	efr = 0.0E0_DP
474
475		! communicate MPI positions
476	<	#ifdef MPI
477	<	call gather(q,qRow,plan_row3)
478	<	call gather(q,qCol,plan_col3)
476	>	#ifdef IS_MPI
477	>	call gather(q,qRow,plan_row3d)
478	>	call gather(q,qCol,plan_col3d)
479		#endif
480
481
#	Line 453 \| Line 487 \| contains
487
488		nlist = 0
489
490	<
490	>
491
492		do i = 1, nrow
493		point(i) = nlist + 1
494	<	#ifdef MPI
494	>	#ifdef IS_MPI
495		ljAtype_i => identPtrListRow(i)%this
496		tag_i = tagRow(i)
497		rxi = qRow(1,i)
#	Line 472 \| Line 506 \| contains
506		#endif
507
508		inner: do j = j_start, ncol
509	<	#ifdef MPI
509	>	#ifdef IS_MPI
510		! Assign identity pointers and tags
511		ljAtype_j => identPtrListColumn(j)%this
512	<	tag_j = tagCol(j)
512	>	tag_j = tagColumn(j)
513		if (newtons_thrd) then
514		if (tag_i <= tag_j) then
515		if (mod(tag_i + tag_j,2) == 0) cycle inner
#	Line 496 \| Line 530 \| contains
530		#endif
531		rijsq = rxijrxij + ryijryij + rzij*rzij
532
533	<	#ifdef MPI
533	>	#ifdef IS_MPI
534		if (rijsq <= rlistsq .AND. &
535		tag_j /= tag_i) then
536		#else
537	+
538		if (rijsq < rlistsq) then
539		#endif
540
#	Line 510 \| Line 545 \| contains
545		endif
546		list(nlist) = j
547
548	<
548	>
549		if (rijsq < rcutsq) then
550
551		r = dsqrt(rijsq)
552
553		call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
554
555	<	#ifdef MPI
555	>	#ifdef IS_MPI
556		eRow(i) = eRow(i) + pot*0.5
557		eCol(i) = eCol(i) + pot*0.5
558		#else
559	<	pe = pe + pot
559	>	pe = pe + pot
560		#endif
561
562	<	efr(1,j) = -rxij
563	<	efr(2,j) = -ryij
564	<	efr(3,j) = -rzij
565	<
566	<	do dim = 1, 3
567	<
568	<
569	<	drdx1 = efr(dim,j) / r
570	<	ftmp = dudr * drdx1
571	<
572	<
573	<	#ifdef MPI
574	<	fCol(dim,j) = fCol(dim,j) - ftmp
575	<	fRow(dim,i) = fRow(dim,i) + ftmp
576	<	#else
577	<
578	<	f(dim,j) = f(dim,j) - ftmp
579	<	f(dim,i) = f(dim,i) + ftmp
580	<
581	<	#endif
582	<	enddo
583	<	endif
584	<	endif
585	<	enddo inner
562	>	drdx = -rxij / r
563	>	drdy = -ryij / r
564	>	drdz = -rzij / r
565	>
566	>	fx = dudr * drdx
567	>	fy = dudr * drdy
568	>	fz = dudr * drdz
569	>
570	>	#ifdef IS_MPI
571	>	fCol(1,j) = fCol(1,j) - fx
572	>	fCol(2,j) = fCol(2,j) - fy
573	>	fCol(3,j) = fCol(3,j) - fz
574	>
575	>	fRow(1,j) = fRow(1,j) + fx
576	>	fRow(2,j) = fRow(2,j) + fy
577	>	fRow(3,j) = fRow(3,j) + fz
578	>	#else
579	>	f(1,j) = f(1,j) - fx
580	>	f(2,j) = f(2,j) - fy
581	>	f(3,j) = f(3,j) - fz
582	>	f(1,i) = f(1,i) + fx
583	>	f(2,i) = f(2,i) + fy
584	>	f(3,i) = f(3,i) + fz
585	>	#endif
586	>
587	>	if (do_stress) then
588	>	tauTemp(1) = tauTemp(1) + fx * rxij
589	>	tauTemp(2) = tauTemp(2) + fx * ryij
590	>	tauTemp(3) = tauTemp(3) + fx * rzij
591	>	tauTemp(4) = tauTemp(4) + fy * rxij
592	>	tauTemp(5) = tauTemp(5) + fy * ryij
593	>	tauTemp(6) = tauTemp(6) + fy * rzij
594	>	tauTemp(7) = tauTemp(7) + fz * rxij
595	>	tauTemp(8) = tauTemp(8) + fz * ryij
596	>	tauTemp(9) = tauTemp(9) + fz * rzij
597	>	endif
598	>	endif
599	>	enddo inner
600		enddo
601
602	<	#ifdef MPI
602	>	#ifdef IS_MPI
603		point(nrow + 1) = nlist + 1
604		#else
605		point(natoms) = nlist + 1
#	Line 564 \| Line 613 \| contains
613		JEND = POINT(i+1) - 1
614		! check thiat molecule i has neighbors
615		if (jbeg .le. jend) then
616	<	#ifdef MPI
616	>	#ifdef IS_MPI
617		ljAtype_i => identPtrListRow(i)%this
618		rxi = qRow(1,i)
619		ryi = qRow(2,i)
#	Line 577 \| Line 626 \| contains
626		#endif
627		do jnab = jbeg, jend
628		j = list(jnab)
629	<	#ifdef MPI
629	>	#ifdef IS_MPI
630		ljAtype_j = identPtrListColumn(j)%this
631	<	rxij = wrap(rxi - q_col(1,j), 1)
632	<	ryij = wrap(ryi - q_col(2,j), 2)
633	<	rzij = wrap(rzi - q_col(3,j), 3)
631	>	rxij = wrap(rxi - qCol(1,j), 1)
632	>	ryij = wrap(ryi - qCol(2,j), 2)
633	>	rzij = wrap(rzi - qCol(3,j), 3)
634		#else
635		ljAtype_j = identPtrList(j)%this
636		rxij = wrap(rxi - q(1,j), 1)
#	Line 595 \| Line 644 \| contains
644		r = dsqrt(rijsq)
645
646		call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
647	<	#ifdef MPI
647	>	#ifdef IS_MPI
648		eRow(i) = eRow(i) + pot*0.5
649		eCol(i) = eCol(i) + pot*0.5
650		#else
651		pe = pe + pot
652		#endif
653	<
654	<
655	<	efr(1,j) = -rxij
656	<	efr(2,j) = -ryij
657	<	efr(3,j) = -rzij
658	<
659	<	do dim = 1, 3
660	<
661	<	drdx1 = efr(dim,j) / r
662	<	ftmp = dudr * drdx1
663	<	#ifdef MPI
664	<	fCol(dim,j) = fCol(dim,j) - ftmp
665	<	fRow(dim,i) = fRow(dim,i) + ftmp
666	<	#else
667	<	f(dim,j) = f(dim,j) - ftmp
668	<	f(dim,i) = f(dim,i) + ftmp
669	<	#endif
670	<	enddo
671	<	endif
672	<	enddo
673	<	endif
674	<	enddo
675	<	endif
653	>
654	>	drdx = -rxij / r
655	>	drdy = -ryij / r
656	>	drdz = -rzij / r
657	>
658	>	fx = dudr * drdx
659	>	fy = dudr * drdy
660	>	fz = dudr * drdz
661	>
662	>	#ifdef IS_MPI
663	>	fCol(1,j) = fCol(1,j) - fx
664	>	fCol(2,j) = fCol(2,j) - fy
665	>	fCol(3,j) = fCol(3,j) - fz
666	>
667	>	fRow(1,j) = fRow(1,j) + fx
668	>	fRow(2,j) = fRow(2,j) + fy
669	>	fRow(3,j) = fRow(3,j) + fz
670	>	#else
671	>	f(1,j) = f(1,j) - fx
672	>	f(2,j) = f(2,j) - fy
673	>	f(3,j) = f(3,j) - fz
674	>	f(1,i) = f(1,i) + fx
675	>	f(2,i) = f(2,i) + fy
676	>	f(3,i) = f(3,i) + fz
677	>	#endif
678	>
679	>	if (do_stress) then
680	>	tauTemp(1) = tauTemp(1) + fx * rxij
681	>	tauTemp(2) = tauTemp(2) + fx * ryij
682	>	tauTemp(3) = tauTemp(3) + fx * rzij
683	>	tauTemp(4) = tauTemp(4) + fy * rxij
684	>	tauTemp(5) = tauTemp(5) + fy * ryij
685	>	tauTemp(6) = tauTemp(6) + fy * rzij
686	>	tauTemp(7) = tauTemp(7) + fz * rxij
687	>	tauTemp(8) = tauTemp(8) + fz * ryij
688	>	tauTemp(9) = tauTemp(9) + fz * rzij
689	>	endif
690	>
691	>
692	>	endif
693	>	enddo
694	>	endif
695	>	enddo
696	>	endif
697	>
698
699
700	<
630	<	#ifdef MPI
700	>	#ifdef IS_MPI
701		!!distribute forces
632	–	call scatter(fRow,f,plan_row3)
702
703	<	call scatter(fCol,f_tmp,plan_col3)
703	>	call scatter(fRow,f,plan_row3d)
704
705	+	call scatter(fCol,fMPITemp,plan_col3d)
706	+
707		do i = 1,nlocal
708	<	f(1:3,i) = f(1:3,i) + f_tmp(1:3,i)
708	>	f(1:3,i) = f(1:3,i) + fMPITemp(1:3,i)
709		end do
710
711
712
713		if (do_pot) then
714	<	! scatter/gather pot_row into the members of my column
715	<	call scatter(e_row,e_tmp,plan_row)
714	>	! scatter/gather pot_row into the members of my column
715	>	call scatter(eRow,eTemp,plan_row)
716
717		! scatter/gather pot_local into all other procs
718		! add resultant to get total pot
719		do i = 1, nlocal
720	<	pe_local = pe_local + e_tmp(i)
720	>	pe_local = pe_local + eTemp(i)
721		enddo
722		if (newtons_thrd) then
723	<	e_tmp = 0.0E0_DP
724	<	call scatter(e_col,e_tmp,plan_col)
723	>	eTemp = 0.0E0_DP
724	>	call scatter(eCol,eTemp,plan_col)
725		do i = 1, nlocal
726	<	pe_local = pe_local + e_tmp(i)
726	>	pe_local = pe_local + eTemp(i)
727		enddo
728		endif
729	+	pe = pe_local
730		endif
731	<	potE = pe_local
731	>
732		#endif
733
734		potE = pe
735
736
737	+	if (do_stress) then
738	+	#ifdef IS_MPI
739	+	mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
740	+	mpi_comm_world,mpi_err)
741	+	#else
742	+	tau = tauTemp
743	+	#endif
744	+	endif
745	+
746	+
747		end subroutine do_lj_ff
748
749		!! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
#	Line 717 \| Line 799 \| contains
799		epsilon = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
800
801
802	<
802	>
803	>
804		call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat)
805
806		r2 = r * r
#	Line 756 \| Line 839 \| contains
839		real(kind = dp) :: param1
840		real(kind = dp) :: param2
841		real(kind = dp ) :: myMixParam
842	+	character(len = getStringLen()) :: thisMixingRule
843		integer, optional :: status
844
845	<
845	>	!! get the mixing rules from the simulation
846	>	thisMixingRule = returnMixingRules()
847		myMixParam = 0.0_dp
848
849		if (present(status)) status = 0
850	<
851	<	select case (thisParam)
852	<
853	<	case ("sigma")
854	<	myMixParam = 0.5_dp * (param1 + param2)
855	<	case ("epsilon")
856	<	myMixParam = sqrt(param1 * param2)
850	>	select case (thisMixingRule)
851	>	case ("standard")
852	>	select case (thisParam)
853	>	case ("sigma")
854	>	myMixParam = 0.5_dp * (param1 + param2)
855	>	case ("epsilon")
856	>	myMixParam = sqrt(param1 * param2)
857	>	case default
858	>	status = -1
859	>	end select
860	>	case("LJglass")
861		case default
862		status = -1
863		end select
775	–
864		end function calcLJMix
865
866

Diff Legend

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

Comparing trunk/mdtools/md_code/lj_FF.F90 (file contents): Revision 252 by chuckv, Tue Jan 28 22:16:55 2003 UTC vs. Revision 281 by chuckv, Mon Feb 24 21:25:15 2003 UTC

Diff Legend

Comparing trunk/mdtools/md_code/lj_FF.F90 (file contents):
Revision 252 by chuckv, Tue Jan 28 22:16:55 2003 UTC vs.
Revision 281 by chuckv, Mon Feb 24 21:25:15 2003 UTC