--- trunk/mdtools/md_code/lj_FF.F90 2003/01/02 21:45:45 222 +++ trunk/mdtools/md_code/lj_FF.F90 2003/01/09 19:40:38 230 @@ -1,42 +1,106 @@ module lj_ff use simulation use definitions, ONLY : dp, ndim +#ifdef IS_MPI + use mpiSimulation +#endif implicit none + PRIVATE +!! Number of lj_atypes in lj_atype_list integer, save :: n_lj_atypes = 0 +!! Starting Size for ljMixed Array + integer, parameter :: ljMixed_blocksize = 10 + type, public :: lj_atype private sequence - integer :: atype_ident = 0 +!! Unique number for place in linked list + integer :: atype_number = 0 +!! Name of atype character(len = 15) :: name = "NULL" +!! Unique indentifier number (ie atomic no, etc) + integer :: atype_ident = 0 +!! Mass of Particle real ( kind = dp ) :: mass = 0.0_dp +!! Lennard-Jones epslon real ( kind = dp ) :: epslon = 0.0_dp +!! Lennard-Jones Sigma real ( kind = dp ) :: sigma = 0.0_dp +!! Lennard-Jones Sigma Squared + real ( kind = dp ) :: sigma2 = 0.0_dp +!! Lennard-Jones Sigma to sixth + real ( kind = dp ) :: sigma6 = 0.0_dp +!! Pointer for linked list creation type (lj_atype), pointer :: next => null() end type lj_atype +!! Pointer type for atype ident array + type :: lj_atypePtr + type (lj_atype), pointer :: this => null() + end type lj_atypePtr + +!! 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() +!! identity pointer list for force loop. + type (lj_atypePtr), dimension(:), allocatable :: identPtrList + +!! Neighbor list and commom arrays + integer, allocatable, dimension(:) :: point + integer, allocatable, dimension(:) :: list + +#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. + real( kind = dp ), allocatable, dimension(:,:) :: qRow + real( kind = dp ), allocatable, dimension(:,:) :: qColumn + + real( kind = dp ), allocatable, dimension(:,:) :: fRow + real( kind = dp ), allocatable, dimension(:,:) :: fColumn +#endif + + + + + + +!! Public methods and data public :: new_lj_atype + public :: do_lj_ff + public :: getLjPot + + + + contains - subroutine new_lj_atype(name,mass,epslon,sigma,status) + subroutine new_lj_atype(name,ident,mass,epslon,sigma,status) character( len = 15 ) :: name real( kind = dp ), intent(in) :: mass real( kind = dp ), intent(in) :: epslon - real( kind = dp ), intent(in) :: sigam + real( kind = dp ), intent(in) :: sigma + integer, intent(in) :: ident integer, intent(out) :: status 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 integer :: alloc_error integer :: atype_counter = 0 + integer :: alloc_size status = 0 + + +! allocate a new atype allocate(this_lj_atype,stat=alloc_error) if (alloc_error /= 0 ) then status = -1 @@ -48,46 +112,580 @@ contains 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%sigma2 +! assume that this atype will be successfully added + this_lj_atype%atype_ident = ident + this_lj_atype%number = n_lj_atypes + 1 +! First time through allocate a array of size ljMixed_blocksize + if(.not. associated(ljMixed)) then + allocate(thisMix(ljMixed_blocksize,ljMixed_blocksize)) + if (alloc_error /= 0 ) then + status = -1 + return + end if + ljMixed => thisMix +! 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) + allocate(thisMix(alloc_size,alloc_size)) + if (alloc_error /= 0 ) then + status = -1 + return + end if +! point oldMix at old ljMixed array + oldMix => ljMixed +! Copy oldMix into new Mixed array + thisMix = oldMix +! Point ljMixed at new array + ljMixed => thisMix +! Free old array so we don't have a memory leak + deallocate(oldMix) + endif + + + + + +! 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 - else ! we need to find the bottom of the list + + else ! we need to find the bottom of the list to insert new atype lj_atype_ptr => lj_atype_list%next + atype_counter = 1 find_end: do if (.not. associated(lj_atype_ptr%next)) then exit find_end end if - lj_atype_ptr => lj_atype_ptr%next +! Set up mixing for new atype and current atype in list + ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma = & + calcLJMix("sigma",this_lj_atype%sigma, & + lj_atype_prt%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 + + 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 + + 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 + end do find_end end if + + + +! 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 - subroutine add_lj_atype() + subroutine init_ljFF() - end subroutine add_lj_atype +#ifdef IS_MPI + - subroutine do_lj_ff(q,f,potE) +#endif + + end subroutine init_ljFF + +!! Takes an ident array and creates an atype pointer list +!! based on those identities + subroutine new_ljatypePtrList(mysize,ident,PtrList,status) + integer, intent(in) :: mysize + integer, intent(in) :: ident + integer, optional :: status + type(lj_atypePtr), dimension(:) :: PtrList + + integer :: thisIdent + integer :: i + integer :: alloc_error + type (lj_atype), pointer :: tmpPtr + + if (present(status)) status = 0 + +! First time through, allocate list + if (.not.(allocated)) then + allocate(PtrList(mysize)) + else +! We want to creat a new ident list so free old list + deallocate(PrtList) + allocate(PtrList(mysize)) + endif + +! Match pointer list + do i = 1, mysize + thisIdent = ident(i) + call getLJatype(thisIdent,tmpPtr) + + if (.not. associated(tmpPtr)) then + status = -1 + return + endif + + PtrList(i)%this => tmpPtr + end do + + end subroutine new_ljatypePtrList + +!! Finds a lj_atype based upon numerical ident +!! 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 :: tmplj_atype_ptr => null() + + if(.not. associated(lj_atype_list)) return + +! Point at head of list. + tmplj_atype_ptr => lj_atype_list + find_ident: do + if (.not.associated(tmplj_atype_ptr)) then + exit find_ident + else if( lj_atype_ptr%atype_ident == ident) + ljAtypePtr => tmplj_atype_ptr + exit find_ident + endif + tmplj_atype_ptr => tmplj_atype_ptr%next + end do find_ident + + end subroutine getLJatype + + +! FORCE routine +!-------------------------------------------------------------> + subroutine do_lj_ff(q,f,potE,do_pot) real ( kind = dp ), dimension(ndim,) :: q real ( kind = dp ), dimension(ndim,nLRparticles) :: f real ( kind = dp ) :: potE + +#ifdef MPI + real( kind = DP ), dimension(3,ncol) :: efr + real( kind = DP ) :: pot_local +#else +! real( kind = DP ), dimension(3,natoms) :: efr +#endif + + real( kind = DP ) :: pe + logical, :: update_nlist + logical :: do_pot + + integer :: i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2 + integer :: nlist + integer :: j_start + 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 + + integer :: nrow + integer :: ncol + + + +#ifndef IS_MPI + nrow = natoms - 1 + ncol = natoms +#else + j_start = 1 +#endif + + call check(update_nlist) + + do_pot = .false. + if (present(pe)) do_pot = .true. + +#ifndef IS_MPI + nloops = nloops + 1 + pot = 0.0E0_DP + f = 0.0E0_DP + e = 0.0E0_DP +#else + f_row = 0.0E0_DP + f_col = 0.0E0_DP + + pot_local = 0.0E0_DP + + e_row = 0.0E0_DP + e_col = 0.0E0_DP + e_tmp = 0.0E0_DP +#endif + efr = 0.0E0_DP + +! communicate MPI positions +#ifdef MPI + call gather(q,q_row,plan_row3) + call gather(q,q_col,plan_col3) +#endif + +#ifndef MPI + +#endif + + if (update_nlist) then + + ! save current configuration, contruct neighbor list, + ! and calculate forces + call save_nlist() + + nlist = 0 + + + + do i = 1, nrow + point(i) = nlist + 1 +#ifdef MPI + tag_i = tag_row(i) + rxi = q_row(1,i) + ryi = q_row(2,i) + rzi = q_row(3,i) +#else + j_start = i + 1 + rxi = q(1,i) + ryi = q(2,i) + rzi = q(3,i) +#endif + inner: do j = j_start, ncol +#ifdef MPI + tag_j = tag_col(j) + if (newtons_thrd) then + if (tag_i <= tag_j) then + if (mod(tag_i + tag_j,2) == 0) cycle inner + else + if (mod(tag_i + tag_j,2) == 1) cycle inner + endif + endif + rxij = wrap(rxi - q_col(1,j), 1) + ryij = wrap(ryi - q_col(2,j), 2) + rzij = wrap(rzi - q_col(3,j), 3) +#else + + rxij = wrap(rxi - q(1,j), 1) + ryij = wrap(ryi - q(2,j), 2) + rzij = wrap(rzi - q(3,j), 3) + +#endif + rijsq = rxij*rxij + ryij*ryij + rzij*rzij +#ifdef MPI + if (rijsq <= rlstsq .AND. & + tag_j /= tag_i) then +#else + if (rijsq < rlstsq) then +#endif + + nlist = nlist + 1 + if (nlist > size(list)) then + call info("FORCE_LJ","error size list smaller then nlist") + write(msg,*) "nlist size(list)", nlist, size(list) + call info("FORCE_LJ",msg) +#ifdef MPI + call mpi_abort(MPI_COMM_WORLD,mpi_err) +#endif + stop + endif + list(nlist) = j + + if (rijsq < rcutsq) then + + r = dsqrt(rijsq) + + call LJ_mix(r,pot,dudr,d2,i,j) + +#ifdef MPI + e_row(i) = e_row(i) + pot*0.5 + e_col(i) = e_col(i) + pot*0.5 +#else + pe = pe + pot +#endif + + efr(1,j) = -rxij + efr(2,j) = -ryij + efr(3,j) = -rzij + do dim = 1, 3 + + drdx1 = efr(dim,j) / r + ftmp = dudr * drdx1 + +#ifdef MPI + f_col(dim,j) = f_col(dim,j) - ftmp + f_row(dim,i) = f_row(dim,i) + ftmp +#else + + f(dim,j) = f(dim,j) - ftmp + f(dim,i) = f(dim,i) + ftmp + +#endif + enddo + endif + endif + enddo inner + enddo + +#ifdef MPI + point(nrow + 1) = nlist + 1 +#else + point(natoms) = nlist + 1 +#endif + + else + + ! use the list to find the neighbors + do i = 1, nrow + JBEG = POINT(i) + JEND = POINT(i+1) - 1 + ! check thiat molecule i has neighbors + if (jbeg .le. jend) then +#ifdef MPI + rxi = q_row(1,i) + ryi = q_row(2,i) + rzi = q_row(3,i) +#else + rxi = q(1,i) + ryi = q(2,i) + rzi = q(3,i) +#endif + do jnab = jbeg, jend + j = list(jnab) +#ifdef MPI + rxij = wrap(rxi - q_col(1,j), 1) + ryij = wrap(ryi - q_col(2,j), 2) + rzij = wrap(rzi - q_col(3,j), 3) +#else + rxij = wrap(rxi - q(1,j), 1) + ryij = wrap(ryi - q(2,j), 2) + rzij = wrap(rzi - q(3,j), 3) +#endif + rijsq = rxij*rxij + ryij*ryij + rzij*rzij + + if (rijsq < rcutsq) then + + r = dsqrt(rijsq) + + call LJ_mix(r,pot,dudr,d2,i,j) +#ifdef MPI + e_row(i) = e_row(i) + pot*0.5 + e_col(i) = e_col(i) + pot*0.5 +#else + if (do_pot) pe = pe + pot +#endif + + + efr(1,j) = -rxij + efr(2,j) = -ryij + efr(3,j) = -rzij + + do dim = 1, 3 + + drdx1 = efr(dim,j) / r + ftmp = dudr * drdx1 +#ifdef MPI + f_col(dim,j) = f_col(dim,j) - ftmp + f_row(dim,i) = f_row(dim,i) + ftmp +#else + f(dim,j) = f(dim,j) - ftmp + f(dim,i) = f(dim,i) + ftmp +#endif + enddo + endif + enddo + endif + enddo + endif + + + +#ifdef MPI + !!distribute forces + call scatter(f_row,f,plan_row3) + + call scatter(f_col,f_tmp,plan_col3) + do i = 1,nlocal + do dim = 1,3 + f(dim,i) = f(dim,i) + f_tmp(dim,i) + end do + end do + + + + if (do_pot) then +! scatter/gather pot_row into the members of my column + call scatter(e_row,e_tmp,plan_row) + + ! 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) + 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) + enddo + endif + endif +#endif + + + + end subroutine do_lj_ff +!! Calculates the potential between two lj particles, optionally returns second +!! derivatives. + subroutine getLjPot(r,pot,dudr,atype1,atype2,d2,status) +! arguments +!! Length of vector between particles + real( kind = dp ), intent(in) :: r +!! Potential Energy + real( kind = dp ), intent(out) :: pot +!! Derivatve wrt postion + real( kind = dp ), intent(out) :: dudr +!! Second Derivative, optional, used mainly for normal mode calculations. + real( kind = dp ), intent(out), optional :: d2 + + type (lj_atype), intent(in) :: atype1 + type (lj_atype), intent(in) :: atype2 + integer, intent(out), optional :: status +! local Variables + real( kind = dp ) :: sigma + real( kind = dp ) :: sigma2 + real( kind = dp ) :: sigma6 + real( kind = dp ) :: epslon + + real( kind = dp ) :: rcut + real( kind = dp ) :: rcut2 + real( kind = dp ) :: rcut6 + real( kind = dp ) :: r2 + real( kind = dp ) :: r6 + + real( kind = dp ) :: t6 + real( kind = dp ) :: t12 + real( kind = dp ) :: tp6 + real( kind = dp ) :: tp12 + real( kind = dp ) :: delta + + logical :: doSec = .false. + + integer :: errorStat + +!! Optional Argument Checking +! Check to see if we need to do second derivatives + + if (present(d2)) doSec = .true. + 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 + + + + call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat) + + r2 = r * r + r6 = r2 * r2 * r2 + + t6 = sigma6/ r6 + t12 = t6 * t6 + + + + tp6 = sigma6 / rcut6 + tp12 = tp6*tp6 + + delta = -4.0E0_DP*epsilon * (tp12 - tp6) + + if (r.le.rcut) then + u = 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 + dudr = 0.0E0_DP + if(doSec) d2 = 0.0E0_DP + endif + + return + + + + end subroutine getLjPot + + +!! Calculates the mixing for sigma or epslon based on character string + function calcLJMix(thisParam,param1,param2,status) result(myMixParam) + character(len=*) :: thisParam + real(kind = dp) :: param1 + real(kind = dp) :: param2 + real(kind = dp ) :: myMixParam + integer, optional :: status + + + myMixParam = 0.0_dp + + if (present(status)) status = 0 + + select case (thisParam) + + case ("sigma") + myMixParam = 0.5_dp * (param1 + param2) + case ("epslon") + myMixParam = sqrt(param1 * param2) + case default + status = -1 + end select + + end function calcLJMix + + + end module lj_ff