--- trunk/mdtools/md_code/lj_FF.F90 2003/01/22 21:45:20 240 +++ trunk/mdtools/md_code/lj_FF.F90 2003/01/24 21:36:52 246 @@ -2,13 +2,13 @@ !! Corresponds to the force field defined in lj_FF.cpp !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: lj_FF.F90,v 1.9 2003-01-22 21:45:20 chuckv Exp $, $Date: 2003-01-22 21:45:20 $, $Name: not supported by cvs2svn $, $Revision: 1.9 $ +!! @version $Id: lj_FF.F90,v 1.10 2003-01-24 21:36:52 chuckv Exp $, $Date: 2003-01-24 21:36:52 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $ module lj_ff use simulation - use definitions, ONLY : dp, ndim + use definitions #ifdef IS_MPI use mpiSimulation #endif @@ -32,7 +32,7 @@ module lj_ff !! Mass of Particle real ( kind = dp ) :: mass = 0.0_dp !! Lennard-Jones epslon - real ( kind = dp ) :: epslon = 0.0_dp + real ( kind = dp ) :: epsilon = 0.0_dp !! Lennard-Jones Sigma real ( kind = dp ) :: sigma = 0.0_dp !! Lennard-Jones Sigma Squared @@ -51,9 +51,9 @@ module lj_ff !! Global list of lj atypes in simulation type (lj_atype), pointer :: lj_atype_list => null() !! LJ mixing array - type (lj_atype), dimension(:,:), allocatable, pointer :: ljMixed =>null() + type (lj_atype), dimension(:,:), pointer :: ljMixed => null() !! identity pointer list for force loop. - type (lj_atypePtr), dimension(:), allocatable :: identPtrList + type (lj_atypePtr), dimension(:), pointer :: identPtrList => null() !! Neighbor list and commom arrays @@ -67,35 +67,47 @@ module lj_ff #ifdef IS_MPI ! Universal routines: All types of force calculations will need these arrays ! Arrays specific to a type of force calculation should be declared in that module. +!! Row position array. real( kind = dp ), allocatable, dimension(:,:) :: qRow +!! Column position array. real( kind = dp ), allocatable, dimension(:,:) :: qColumn +!! Row Force array. real( kind = dp ), allocatable, dimension(:,:) :: fRow +!! Column Force Array. real( kind = dp ), allocatable, dimension(:,:) :: fColumn - type (lj_atypePtr), dimension(:), allocatable :: identPtrListRow - type (lj_atypePtr), dimension(:), allocatable :: identPtrListColumn -#endif +!! Row potential energy array + real( kind = dp ), allocatable, dimension(:) :: eRow +!! Column potential energy array + real( kind = dp ), allocatable, dimension(:) :: eColumn +!! Row lj_atype pointer list + type (lj_atypePtr), dimension(:), pointer :: identPtrListRow => null() +!! Column lj_atype pointer list + type (lj_atypePtr), dimension(:), pointer :: identPtrListColumn => null() +#endif - logical :: isljFFinit = .false. +!! Logical has lj force field module been initialized? + logical, save :: isljFFinit = .false. + !! Public methods and data public :: new_lj_atype public :: do_lj_ff public :: getLjPot - + public :: init_ljFF contains - subroutine new_lj_atype(ident,mass,epslon,sigma,status) + subroutine new_lj_atype(ident,mass,epsilon,sigma,status) real( kind = dp ), intent(in) :: mass - real( kind = dp ), intent(in) :: epslon + real( kind = dp ), intent(in) :: epsilon real( kind = dp ), intent(in) :: sigma integer, intent(in) :: ident integer, intent(out) :: status @@ -103,8 +115,8 @@ contains type (lj_atype), pointer :: this_lj_atype type (lj_atype), pointer :: lj_atype_ptr - type (lj_atype), allocatable, dimension(:,:), pointer :: thisMix - type (lj_atype), allocatable, dimension(:,:), pointer :: oldMix + type (lj_atype), dimension(:,:), pointer :: thisMix + type (lj_atype), dimension(:,:), pointer :: oldMix integer :: alloc_error integer :: atype_counter = 0 integer :: alloc_size @@ -121,15 +133,15 @@ contains end if ! assign our new lj_atype information - this_lj_atype%mass = mass - this_lj_atype%epslon = epslon - this_lj_atype%sigma = sigma - this_lj_atype%sigma2 = sigma * sigma - this_lj_atype%sigma6 = this_lj_atype%sigma2 * this_lj_atype%sigma2 & + this_lj_atype%mass = mass + this_lj_atype%epsilon = epsilon + this_lj_atype%sigma = sigma + this_lj_atype%sigma2 = sigma * sigma + this_lj_atype%sigma6 = this_lj_atype%sigma2 * this_lj_atype%sigma2 & * this_lj_atype%sigma2 ! assume that this atype will be successfully added this_lj_atype%atype_ident = ident - this_lj_atype%number = n_lj_atypes + 1 + this_lj_atype%atype_number = n_lj_atypes + 1 ! First time through allocate a array of size ljMixed_blocksize @@ -143,7 +155,7 @@ contains ! If we have outgrown ljMixed_blocksize, allocate a new matrix twice the size and ! point ljMix at the new matrix. else if( (n_lj_atypes + 1) > size(ljMixed)) then - alloc_size = 2*size(ljMix) + alloc_size = 2*size(ljMixed) allocate(thisMix(alloc_size,alloc_size)) if (alloc_error /= 0 ) then status = -1 @@ -166,64 +178,78 @@ contains ! Find bottom of atype master list ! if lj_atype_list is null then we are at the top of the list. if (.not. associated(lj_atype_list)) then - lj_atype_ptr => this_lj_atype - atype_counter = 1 + write(*,*) "Associating lists for first time" + + + this_lj_atype%atype_number = 1 + + ljMixed(n_lj_atypes,n_lj_atypes)%sigma = this_lj_atype%sigma + + ljMixed(n_lj_atypes,n_lj_atypes)%sigma2 = ljMixed(n_lj_atypes,n_lj_atypes)%sigma & + * ljMixed(n_lj_atypes,n_lj_atypes)%sigma + + ljMixed(n_lj_atypes,n_lj_atypes)%sigma6 = ljMixed(n_lj_atypes,n_lj_atypes)%sigma2 & + * ljMixed(n_lj_atypes,n_lj_atypes)%sigma2 & + * ljMixed(n_lj_atypes,n_lj_atypes)%sigma2 + + ljMixed(n_lj_atypes,n_lj_atypes)%epsilon = this_lj_atype%epsilon + + lj_atype_list => this_lj_atype + + write(*,*) "lj_atype_list first time through -- ident", lj_atype_list%atype_ident else ! we need to find the bottom of the list to insert new atype - lj_atype_ptr => lj_atype_list%next + write(*,*) "Second time through" + lj_atype_ptr => lj_atype_list + write(*,*) "lj_atype_ptr to next" atype_counter = 1 find_end: do - if (.not. associated(lj_atype_ptr%next)) then - exit find_end - end if + if (.not. associated(lj_atype_ptr)) exit find_end + + write(*,*) "Inside find_end-this ident", lj_atype_ptr%atype_ident ! Set up mixing for new atype and current atype in list - ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma = & + ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma = & calcLJMix("sigma",this_lj_atype%sigma, & - lj_atype_prt%sigma) + lj_atype_ptr%sigma) + + ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 = & + ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma & + * ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma - ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 = & - ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma & - * ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma + ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma6 = & + ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 & + * ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 & + * ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 - ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma6 = & - ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 & - * ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 & - * ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 + ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%epsilon = & + calcLJMix("epsilon",this_lj_atype%epsilon, & + lj_atype_ptr%epsilon) - ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%epslon = & - calcLJMix("epslon",this_lj_atype%epslon, & - lj_atype_prt%epslon) - ! Advance to next pointer - lj_atype_ptr => lj_atype_ptr%next - atype_counter = atype_counter + 1 - + lj_atype_ptr => lj_atype_ptr%next + atype_counter = atype_counter + 1 end do find_end + + lj_atype_ptr => this_lj_atype + write(*,*) "just added: ", lj_atype_ptr%atype_ident end if - - + write(*,*) "lj_atype_list now contains.." + lj_atype_ptr => lj_atype_list + do while (associated(lj_atype_ptr)) + write(*,*) "lj_atype_list contains ident", lj_atype_ptr%atype_ident + lj_atype_ptr => lj_atype_ptr%next + end do + ! Insert new atype at end of list - lj_atype_ptr => this_lj_atype + ! Increment number of atypes n_lj_atypes = n_lj_atypes + 1 ! Set up self mixing - ljMix(n_lj_atypes,n_lj_atypes)%sigma = this_lj_atype%sigma - - ljMix(n_lj_atypes,n_lj_atypes)%sigma2 = ljMix(n_lj_atypes,n_lj_atypes)%sigma & - * ljMix(n_lj_atypes,n_lj_atypes)%sigma - - ljMix(n_lj_atypes,n_lj_atypes)%sigma6 = ljMix(n_lj_atypes,n_lj_atypes)%sigma2 & - * ljMix(n_lj_atypes,n_lj_atypes)%sigma2 & - * ljMix(n_lj_atypes,n_lj_atypes)%sigma2 - - ljMix(n_lj_atypes,n_lj_atypes)%epslon = this_lj_atype%epslon - - end subroutine new_lj_atype @@ -236,7 +262,10 @@ contains !! Result status, success = 0, error = -1 integer, intent(out) :: Status + integer :: alloc_stat + integer :: thisStat + integer :: i #ifdef IS_MPI integer, allocatable, dimension(:) :: identRow integer, allocatable, dimension(:) :: identCol @@ -246,6 +275,7 @@ contains #endif status = 0 + write(*,*) "Initializing ljFF" !! if were're not in MPI, we just update ljatypePtrList #ifndef IS_MPI call new_ljatypePtrList(nComponents,ident,identPtrList,thisStat) @@ -253,6 +283,7 @@ contains status = -1 return endif + write(*,*) "component pointer list initialized" !! Allocate pointer lists allocate(point(nComponents),stat=alloc_stat) if (alloc_stat /=0) then @@ -371,6 +402,53 @@ contains return endif endif + + if (.not.allocated(eRow)) then + allocate(eRow(nrow),stat=alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif + else + deallocate(eRow) + allocate(eRow(nrow),stat=alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif + endif + + if (.not.allocated(eCol)) then + allocate(eCol(ncol),stat=alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif + else + deallocate(eCol) + allocate(eCol(ncol),stat=alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif + endif + + if (.not.allocated(eTemp)) then + allocate(eTemp(nComponents),stat=alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif + else + deallocate(eTemp) + allocate(eTemp(nComponets),stat=alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif + endif + + !! Allocate neighbor lists for mpi simulations. if (.not. allocated(point)) then allocate(point(nrow),stat=alloc_stat) @@ -402,10 +480,13 @@ contains #endif + + do i = 1, nComponents + write(*,*) "Ident pointer list ident: ", identPtrList(i)%this%atype_ident + end do isljFFinit = .true. - end subroutine init_ljFF @@ -416,9 +497,9 @@ contains !! based on those identities subroutine new_ljatypePtrList(mysize,ident,PtrList,status) integer, intent(in) :: mysize - integer, intent(in) :: ident + integer,dimension(:), intent(in) :: ident integer, optional :: status - type(lj_atypePtr), dimension(:) :: PtrList + type(lj_atypePtr), dimension(:), pointer :: PtrList integer :: thisIdent integer :: i @@ -427,21 +508,26 @@ contains if (present(status)) status = 0 + write(*,*) "Creating new ljatypePtrList...." ! First time through, allocate list - if (.not.(allocated)) then + if (.not.associated(PtrList)) then + write(*,*) "allocating pointer list...." allocate(PtrList(mysize)) else ! We want to creat a new ident list so free old list - deallocate(PrtList) + deallocate(PtrList) allocate(PtrList(mysize)) endif ! Match pointer list + write(*,*) "Doing ident search" do i = 1, mysize thisIdent = ident(i) + write(*,*) "Calling getLJatype for ident ", thisIdent call getLJatype(thisIdent,tmpPtr) if (.not. associated(tmpPtr)) then + write(*,*) "tmpptr was not associated" status = -1 return endif @@ -455,32 +541,41 @@ contains !! returns a null pointer if error subroutine getLJatype(ident,ljAtypePtr) integer, intent(in) :: ident - type (lj_atype), intent(out),pointer :: ljAtypePtr => null() - + type (lj_atype), pointer :: ljAtypePtr type (lj_atype), pointer :: tmplj_atype_ptr => null() + write(*,*) "Finding ljAtype for ident ", ident + ljAtypePtr => null() + if(.not. associated(lj_atype_list)) return ! Point at head of list. + write(*,*) "Searching at head of list" tmplj_atype_ptr => lj_atype_list find_ident: do if (.not.associated(tmplj_atype_ptr)) then + write(*,*) "Find_ident, pointer not associated" exit find_ident - else if( lj_atype_ptr%atype_ident == ident) + else if( tmplj_atype_ptr%atype_ident == ident) then ljAtypePtr => tmplj_atype_ptr exit find_ident endif tmplj_atype_ptr => tmplj_atype_ptr%next end do find_ident + + + end subroutine getLJatype !! FORCE routine Calculates Lennard Jones forces. !-------------------------------------------------------------> subroutine do_lj_ff(q,f,potE,do_pot) - real ( kind = dp ), dimension(ndim,) :: q - real ( kind = dp ), dimension(ndim,nLRparticles) :: f +!! Position array provided by C, dimensioned by getNlocal + real ( kind = dp ), dimension(3,getNlocal()) :: q +!! Force array provided by C, dimensioned by getNlocal + real ( kind = dp ), dimension(3,getNlocal()) :: f real ( kind = dp ) :: potE logical ( kind = 2) :: do_pot @@ -488,14 +583,14 @@ contains type(lj_atype), pointer :: ljAtype_j #ifdef MPI - real( kind = DP ), dimension(3,ncol) :: efr + real( kind = DP ), dimension(3,getNcol()) :: efr real( kind = DP ) :: pot_local #else -! real( kind = DP ), dimension(3,natoms) :: efr + real( kind = DP ), dimension(3,getNlocal()) :: efr #endif real( kind = DP ) :: pe - logical, :: update_nlist + logical :: update_nlist integer :: i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2 @@ -504,47 +599,66 @@ contains integer :: tag_i,tag_j real( kind = DP ) :: r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp real( kind = DP ) :: rxi, ryi, rzi, rxij, ryij, rzij, rijsq - real( kind = DP ) :: rlistsq, rcutsq + real( kind = DP ) :: rlistsq, rcutsq,rlist,rcut +! a rig that need to be fixed. +#ifdef IS_MPI + logical :: newtons_thrd = .true. + real( kind = dp ) :: pe_local + integer :: nlocal +#endif integer :: nrow integer :: ncol + integer :: natoms +! Make sure we are properly initialized. if (.not. isljFFInit) then write(default_error,*) "ERROR: lj_FF has not been properly initialized" + return + endif +#ifdef IS_MPI + if (.not. isMPISimSet()) then + write(default_error,*) "ERROR: mpiSimulation has not been properly initialized" return endif +#endif +!! initialize local variables + natoms = getNlocal() + call getRcut(rcut,rcut2=rcutsq) + call getRlist(rlist,rlistsq) #ifndef IS_MPI nrow = natoms - 1 ncol = natoms #else nrow = getNrow(plan_row) ncol = getNcol(plan_col) + nlocal = natoms j_start = 1 #endif +!! See if we need to update neighbor lists + call check(q,update_nlist) - call check(update_nlist) - !--------------WARNING........................... ! Zero variables, NOTE:::: Forces are zeroed in C ! Zeroing them here could delete previously computed ! Forces. !------------------------------------------------ #ifndef IS_MPI - nloops = nloops + 1 - pot = 0.0E0_DP - e = 0.0E0_DP +! nloops = nloops + 1 + pe = 0.0E0_DP + #else f_row = 0.0E0_DP f_col = 0.0E0_DP - pot_local = 0.0E0_DP + pe_local = 0.0E0_DP - e_row = 0.0E0_DP - e_col = 0.0E0_DP - e_tmp = 0.0E0_DP + eRow = 0.0E0_DP + eCol = 0.0E0_DP + eTemp = 0.0E0_DP #endif efr = 0.0E0_DP @@ -554,15 +668,12 @@ contains call gather(q,qCol,plan_col3) #endif -#ifndef MPI -#endif - if (update_nlist) then ! save current configuration, contruct neighbor list, ! and calculate forces - call save_nlist() + call save_nlist(q) nlist = 0 @@ -610,15 +721,15 @@ contains rijsq = rxij*rxij + ryij*ryij + rzij*rzij #ifdef MPI - if (rijsq <= rlstsq .AND. & + if (rijsq <= rlistsq .AND. & tag_j /= tag_i) then #else - if (rijsq < rlstsq) then + if (rijsq < rlistsq) then #endif nlist = nlist + 1 if (nlist > size(list)) then -#warning "Change how nlist size is done" +!! "Change how nlist size is done" write(DEFAULT_ERROR,*) "ERROR: nlist > list size" endif list(nlist) = j @@ -631,8 +742,8 @@ contains call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j) #ifdef MPI - e_row(i) = e_row(i) + pot*0.5 - e_col(i) = e_col(i) + pot*0.5 + eRow(i) = eRow(i) + pot*0.5 + eCol(i) = eCol(i) + pot*0.5 #else pe = pe + pot #endif @@ -709,10 +820,10 @@ contains call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j) #ifdef MPI - e_row(i) = e_row(i) + pot*0.5 - e_col(i) = e_col(i) + pot*0.5 + eRow(i) = eRow(i) + pot*0.5 + eCol(i) = eCol(i) + pot*0.5 #else - if (do_pot) pe = pe + pot + pe = pe + pot #endif @@ -745,10 +856,9 @@ contains call scatter(fRow,f,plan_row3) call scatter(fCol,f_tmp,plan_col3) + do i = 1,nlocal - do dim = 1,3 - f(dim,i) = f(dim,i) + f_tmp(dim,i) - end do + f(1:3,i) = f(1:3,i) + f_tmp(1:3,i) end do @@ -760,21 +870,22 @@ contains ! scatter/gather pot_local into all other procs ! add resultant to get total pot do i = 1, nlocal - pot_local = pot_local + e_tmp(i) + pe_local = pe_local + e_tmp(i) enddo if (newtons_thrd) then e_tmp = 0.0E0_DP call scatter(e_col,e_tmp,plan_col) do i = 1, nlocal - pot_local = pot_local + e_tmp(i) + pe_local = pe_local + e_tmp(i) enddo endif endif + potE = pe_local #endif + potE = pe - end subroutine do_lj_ff !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second @@ -790,8 +901,8 @@ contains !! Second Derivative, optional, used mainly for normal mode calculations. real( kind = dp ), intent(out), optional :: d2 - type (lj_atype), intent(in), pointer :: atype1 - type (lj_atype), intent(in), pointer :: atype2 + type (lj_atype), pointer :: atype1 + type (lj_atype), pointer :: atype2 integer, intent(out), optional :: status @@ -799,7 +910,7 @@ contains real( kind = dp ) :: sigma real( kind = dp ) :: sigma2 real( kind = dp ) :: sigma6 - real( kind = dp ) :: epslon + real( kind = dp ) :: epsilon real( kind = dp ) :: rcut real( kind = dp ) :: rcut2 @@ -824,10 +935,10 @@ contains if (present(status)) status = 0 ! Look up the correct parameters in the mixing matrix - sigma = ljMixed(atype1%atype_ident,atype2_atype_ident)%sigma - sigma2 = ljMixed(atype1%atype_ident,atype2_atype_ident)%sigma2 - sigma6 = ljMixed(atype1%atype_ident,atype2_atype_ident)%sigma6 - epslon = ljMixed(atype1%atype_ident,atype2_atype_ident)%epslon + sigma = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma + sigma2 = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma2 + sigma6 = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma6 + epsilon = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon @@ -847,11 +958,11 @@ contains delta = -4.0E0_DP*epsilon * (tp12 - tp6) if (r.le.rcut) then - u = 4.0E0_DP * epsilon * (t12 - t6) + delta + pot = 4.0E0_DP * epsilon * (t12 - t6) + delta dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / r if(doSec) d2 = 24.0E0_DP * epsilon * (26.0E0_DP*t12 - 7.0E0_DP*t6)/r/r else - u = 0.0E0_DP + pot = 0.0E0_DP dudr = 0.0E0_DP if(doSec) d2 = 0.0E0_DP endif @@ -880,7 +991,7 @@ contains case ("sigma") myMixParam = 0.5_dp * (param1 + param2) - case ("epslon") + case ("epsilon") myMixParam = sqrt(param1 * param2) case default status = -1