--- trunk/mdtools/md_code/lj_FF.F90 2003/01/10 21:56:10 235 +++ trunk/mdtools/md_code/lj_FF.F90 2003/01/27 18:28:11 247 @@ -1,6 +1,15 @@ +!! Calculates Long Range forces Lennard-Jones interactions. +!! 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.11 2003-01-27 18:28:11 chuckv Exp $, $Date: 2003-01-27 18:28:11 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $ + + + module lj_ff use simulation - use definitions, ONLY : dp, ndim + use definitions + use generic_atypes #ifdef IS_MPI use mpiSimulation #endif @@ -10,298 +19,357 @@ module lj_ff !! Number of lj_atypes in lj_atype_list integer, save :: n_lj_atypes = 0 -!! Starting Size for ljMixed Array - integer, parameter :: ljMixed_blocksize = 10 +!! Global list of lj atypes in simulation + type (lj_atype), pointer :: ljListHead => null() + type (lj_atype), pointer :: ljListTail => null() - type, public :: lj_atype - private - sequence -!! Unique number for place in linked list - integer :: atype_number = 0 -!! 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, public :: 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 + type (lj_atype), dimension(:,:), pointer :: ljMixed => null() !! Neighbor list and commom arrays +!! This should eventally get moved to a neighbor list type integer, allocatable, dimension(:) :: point integer, allocatable, dimension(:) :: list + integer, parameter :: listMultiplier = 80 + integer :: nListAllocs = 0 + integer, parameter :: maxListAllocs = 5 -#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 +!! Atype identity pointer lists +#ifdef IS_MPI +!! Row lj_atype pointer list + type (lj_atypePtr), dimension(:), pointer :: identPtrListRow => null() +!! Column lj_atype pointer list + type (lj_atypePtr), dimension(:), pointer :: identPtrListColumn => null() +#else + type( lj_identPtrList ), dimension(:), pointer :: identPtrList => null() #endif +!! 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) +!! Adds a new lj_atype to the list. + 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 - 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), pointer :: newLJ_atype integer :: alloc_error integer :: atype_counter = 0 integer :: alloc_size - + integer :: err_stat status = 0 ! allocate a new atype - allocate(this_lj_atype,stat=alloc_error) + allocate(newLJ_atype,stat=alloc_error) if (alloc_error /= 0 ) then status = -1 return 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%sigma2 + newLJ_atype%mass = mass + newLJ_atype%epsilon = epsilon + newLJ_atype%sigma = sigma + newLJ_atype%sigma2 = sigma * sigma + newLJ_atype%sigma6 = newLJ_atype%sigma2 * newLJ_atype%sigma2 & + * newLJ_atype%sigma2 ! assume that this atype will be successfully added - this_lj_atype%atype_ident = ident - this_lj_atype%number = n_lj_atypes + 1 + newLJ_atype%atype_ident = ident + newLJ_atype%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) + call add_atype(newLJ_atype,ljListHead,ljListTail,err_stat) + if (err_stat /= 0 ) then + status = -1 + return endif + n_lj_atypes = n_lj_atypes + 1 + end subroutine new_lj_atype -! 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 + subroutine init_ljFF(nComponents,ident, status) + !! Number of components in ident array + integer, intent(inout) :: nComponents +!! Array of identities nComponents long corresponding to +!! ljatype ident. + integer, dimension(nComponents),intent(inout) :: ident +!! Result status, success = 0, error = -1 + integer, intent(out) :: Status - 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 -! 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) + integer :: alloc_stat - 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 + integer :: thisStat + integer :: i +#ifdef IS_MPI + integer, allocatable, dimension(:) :: identRow + integer, allocatable, dimension(:) :: identCol + integer :: nrow + integer :: ncol + integer :: alloc_stat +#endif + status = 0 +!! if were're not in MPI, we just update ljatypePtrList +#ifndef IS_MPI + call create_IdentPtrlst(ident,ljListHead,identPtrList,thisStat) + if ( thisStat /= 0 ) then + status = -1 + return + endif - 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 +!! Allocate pointer lists + allocate(point(nComponents),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif - 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 - - - + allocate(list(nComponents*listMultiplier),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif -! Insert new atype at end of list - lj_atype_ptr => this_lj_atype -! Increment number of atypes +! if were're in MPI, we also have to worry about row and col lists +#else +! We can only set up forces if mpiSimulation has been setup. + if (.not. isMPISimSet()) then + status = -1 + return + endif + nrow = getNrow() + ncol = getNcol() +!! Allocate temperary arrays to hold gather information + allocate(identRow(nrow),stat=alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif - n_lj_atypes = n_lj_atypes + 1 + allocate(identCol(ncol),stat=alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif +!! Gather idents into row and column idents + call gather(ident,identRow,plan_row) + call gather(ident,identCol,plan_col) -! Set up self mixing +!! Create row and col pointer lists + call create_IdentPtrlst(identRow,ljListTail,identPtrListRow,thisStat) + if (thisStat /= 0 ) then + status = -1 + return + endif - ljMix(n_lj_atypes,n_lj_atypes)%sigma = this_lj_atype%sigma + call create_IdentPtrlst(identCol,ljListTail,identPtrListCol,thisStat) + if (thisStat /= 0 ) then + status = -1 + return + endif - ljMix(n_lj_atypes,n_lj_atypes)%sigma2 = ljMix(n_lj_atypes,n_lj_atypes)%sigma & - * ljMix(n_lj_atypes,n_lj_atypes)%sigma +!! free temporary ident arrays + deallocate(identCol) + deallocate(identRow) - 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 +!! Allocate neighbor lists for mpi simulations. + if (.not. allocated(point)) then + allocate(point(nrow),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif + allocate(list(ncol*listMultiplier),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif + else + deallocate(list) + deallocate(point) + allocate(point(nrow),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif - end subroutine new_lj_atype + allocate(list(ncol*listMultiplier),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif + endif +#endif + + call createMixingList(thisStat) + if (thisStat /= 0) then + status = -1 + return + endif + isljFFinit = .true. - subroutine init_ljFF() + end subroutine init_ljFF -#ifdef IS_MPI - -#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 + subroutine createMixingList(status) + integer :: listSize + integer :: status integer :: i - integer :: alloc_error - type (lj_atype), pointer :: tmpPtr + integer :: j - if (present(status)) status = 0 + integer :: outerCounter = 0 + integer :: innerCounter = 0 + type (lj_atype), pointer :: tmpPtrOuter => null() + type (lj_atype), pointer :: tmpPtrInner => null() + 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 + listSize = getListLen(ljListHead) + if (listSize == 0) then + status = -1 + return + end if + -! Match pointer list - do i = 1, mysize - thisIdent = ident(i) - call getLJatype(thisIdent,tmpPtr) - - if (.not. associated(tmpPtr)) then - status = -1 - return + if (.not. associated(ljMixed)) then + allocate(ljMixed(listSize,listSize)) + else + status = -1 + return + end if + + + write(*,*) "Creating mixing list" + tmpPtrOuter => ljListHead + tmpPtrInner => tmpPtrOuter%next + do while (associated(tmpPtrOuter)) + outerCounter = outerCounter + 1 + write(*,*) "Outer index is: ", outerCounter + write(*,*) "For atype: ", tmpPtrOuter%atype_ident +! do self mixing rule + ljMixed(outerCounter,outerCounter)%sigma = tmpPtrOuter%sigma + + ljMixed(outerCounter,outerCounter)%sigma2 = ljMixed(outerCounter,outerCounter)%sigma & + * ljMixed(outerCounter,outerCounter)%sigma + + ljMixed(outerCounter,outerCounter)%sigma6 = ljMixed(outerCounter,outerCounter)%sigma2 & + * ljMixed(outerCounter,outerCounter)%sigma2 & + * ljMixed(outerCounter,outerCounter)%sigma2 + + ljMixed(outerCounter,outerCounter)%epsilon = tmpPtrOuter%epsilon + + innerCounter = outerCounter + 1 + do while (associated(tmpPtrInner)) + + write(*,*) "Entered inner loop for: ", innerCounter + ljMixed(outerCounter,innerCounter)%sigma = & + calcLJMix("sigma",tmpPtrOuter%sigma, & + tmpPtrInner%sigma) + + ljMixed(outerCounter,innerCounter)%sigma2 = & + ljMixed(outerCounter,innerCounter)%sigma & + * ljMixed(outerCounter,innerCounter)%sigma + + ljMixed(outerCounter,innerCounter)%sigma6 = & + ljMixed(outerCounter,innerCounter)%sigma2 & + * ljMixed(outerCounter,innerCounter)%sigma2 & + * ljMixed(outerCounter,innerCounter)%sigma2 + + ljMixed(outerCounter,innerCounter)%epsilon = & + calcLJMix("epsilon",tmpPtrOuter%epsilon, & + tmpPtrInner%epsilon) + ljMixed(innerCounter,outerCounter)%sigma = ljMixed(outerCounter,innerCounter)%sigma + ljMixed(innerCounter,outerCounter)%sigma2 = ljMixed(outerCounter,innerCounter)%sigma2 + ljMixed(innerCounter,outerCounter)%sigma6 = ljMixed(outerCounter,innerCounter)%sigma6 + ljMixed(innerCounter,outerCounter)%epsilon = ljMixed(outerCounter,innerCounter)%epsilon + + + tmpPtrInner => tmpPtrInner%next + innerCounter = innerCounter + 1 + end do +! advance pointers + tmpPtrOuter => tmpPtrOuter%next + if (associated(tmpPtrOuter)) then + tmpPtrInner => tmpPtrOuter%next endif - PtrList(i)%this => tmpPtr end do - end subroutine new_ljatypePtrList + end subroutine createMixingList -!! 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 + + +!! 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 + + type(lj_atype), pointer :: ljAtype_i + 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 +!! Local arrays needed for MPI +#ifdef IS_MPI + real(kind = dp), dimension(3:getNrow()) :: qRow + real(kind = dp), dimension(3:getNcol()) :: qCol + + real(kind = dp), dimension(3:getNrow()) :: fRow + real(kind = dp), dimension(3:getNcol()) :: fCol + + real(kind = dp), dimension(getNrow()) :: eRow + real(kind = dp), dimension(getNcol()) :: eCol + + real(kind = dp), dimension(getNlocal()) :: eTemp +#endif + + + real( kind = DP ) :: pe - logical, :: update_nlist + logical :: update_nlist integer :: i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2 @@ -310,53 +378,81 @@ 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,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 - call check(update_nlist) + +!! See if we need to update neighbor lists + call check(q,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 - f = 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 ! communicate MPI positions #ifdef MPI - call gather(q,q_row,plan_row3) - call gather(q,q_col,plan_col3) + call gather(q,qRow,plan_row3) + 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 @@ -365,11 +461,13 @@ contains 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) + ljAtype_i => identPtrListRow(i)%this + tag_i = tagRow(i) + rxi = qRow(1,i) + ryi = qRow(2,i) + rzi = qRow(3,i) #else + ljAtype_i => identPtrList(i)%this j_start = i + 1 rxi = q(1,i) ryi = q(2,i) @@ -378,7 +476,9 @@ contains inner: do j = j_start, ncol #ifdef MPI - tag_j = tag_col(j) +! Assign identity pointers and tags + ljAtype_j => identPtrListColumn(j)%this + tag_j = tagCol(j) if (newtons_thrd) then if (tag_i <= tag_j) then if (mod(tag_i + tag_j,2) == 0) cycle inner @@ -387,11 +487,11 @@ contains 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) + rxij = wrap(rxi - qCol(1,j), 1) + ryij = wrap(ryi - qCol(2,j), 2) + rzij = wrap(rzi - qCol(3,j), 3) #else - + ljAtype_j => identPtrList(j)%this rxij = wrap(rxi - q(1,j), 1) ryij = wrap(ryi - q(2,j), 2) rzij = wrap(rzi - q(3,j), 3) @@ -400,21 +500,16 @@ 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 - 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 +!! "Change how nlist size is done" + write(DEFAULT_ERROR,*) "ERROR: nlist > list size" endif list(nlist) = j @@ -423,11 +518,11 @@ contains r = dsqrt(rijsq) - call LJ_mix(r,pot,dudr,d2,i,j) + 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 @@ -444,8 +539,8 @@ contains #ifdef MPI - f_col(dim,j) = f_col(dim,j) - ftmp - f_row(dim,i) = f_row(dim,i) + ftmp + fCol(dim,j) = fCol(dim,j) - ftmp + fRow(dim,i) = fRow(dim,i) + ftmp #else f(dim,j) = f(dim,j) - ftmp @@ -473,10 +568,12 @@ contains ! 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) + ljAtype_i => identPtrListRow(i)%this + rxi = qRow(1,i) + ryi = qRow(2,i) + rzi = qRow(3,i) #else + ljAtype_i => identPtrList(i)%this rxi = q(1,i) ryi = q(2,i) rzi = q(3,i) @@ -484,13 +581,15 @@ contains 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) + ljAtype_j = identPtrListColumn(j)%this + 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) + ljAtype_j = identPtrList(j)%this + 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 @@ -498,12 +597,12 @@ contains r = dsqrt(rijsq) - call LJ_mix(r,pot,dudr,d2,i,j) + 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 @@ -516,8 +615,8 @@ contains 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 + fCol(dim,j) = fCol(dim,j) - ftmp + fRow(dim,i) = fRow(dim,i) + ftmp #else f(dim,j) = f(dim,j) - ftmp f(dim,i) = f(dim,i) + ftmp @@ -533,13 +632,12 @@ contains #ifdef MPI !!distribute forces - call scatter(f_row,f,plan_row3) + call scatter(fRow,f,plan_row3) - call scatter(f_col,f_tmp,plan_col3) + 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 @@ -551,24 +649,25 @@ 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, optionally returns second +!! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second !! derivatives. subroutine getLjPot(r,pot,dudr,atype1,atype2,d2,status) ! arguments @@ -581,8 +680,8 @@ contains !! 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 + type (lj_atype), pointer :: atype1 + type (lj_atype), pointer :: atype2 integer, intent(out), optional :: status @@ -590,7 +689,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 @@ -615,10 +714,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 @@ -638,11 +737,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 @@ -654,7 +753,7 @@ contains end subroutine getLjPot -!! Calculates the mixing for sigma or epslon based on character string +!! Calculates the mixing for sigma or epslon based on character string for initialzition of mixing array. function calcLJMix(thisParam,param1,param2,status) result(myMixParam) character(len=*) :: thisParam real(kind = dp) :: param1 @@ -671,7 +770,7 @@ contains case ("sigma") myMixParam = 0.5_dp * (param1 + param2) - case ("epslon") + case ("epsilon") myMixParam = sqrt(param1 * param2) case default status = -1