--- trunk/mdtools/md_code/lj_FF.F90 2003/01/20 22:36:12 239 +++ trunk/mdtools/md_code/lj_FF.F90 2003/01/30 20:03:37 254 @@ -2,13 +2,14 @@ !! 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.8 2003-01-20 22:36:12 chuckv Exp $, $Date: 2003-01-20 22:36:12 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $ +!! @version $Id: lj_FF.F90,v 1.14 2003-01-30 20:03:36 chuckv Exp $, $Date: 2003-01-30 20:03:36 $, $Name: not supported by cvs2svn $, $Revision: 1.14 $ module lj_ff use simulation - use definitions, ONLY : dp, ndim + use definitions + use generic_atypes #ifdef IS_MPI use mpiSimulation #endif @@ -18,42 +19,13 @@ 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() -!! Basic atom type for a Lennard-Jones Atom. - 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 @@ -64,166 +36,78 @@ module lj_ff 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 + logical, save :: firstTime = .True. - real( kind = dp ), allocatable, dimension(:,:) :: fRow - real( kind = dp ), allocatable, dimension(:,:) :: fColumn - - type (lj_atypePtr), dimension(:), allocatable :: identPtrListRow - type (lj_atypePtr), dimension(:), allocatable :: identPtrListColumn +!! Atype identity pointer lists +#ifdef IS_MPI +!! Row lj_atype pointer list + type (lj_identPtrList), dimension(:), pointer :: identPtrListRow => null() +!! Column lj_atype pointer list + type (lj_identPtrList), 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. - logical :: 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 - - - - - -! 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 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) - - 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 @@ -236,23 +120,29 @@ 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 integer :: nrow integer :: ncol - integer :: alloc_stat #endif status = 0 + + + !! if were're not in MPI, we just update ljatypePtrList #ifndef IS_MPI - call new_ljatypePtrList(nComponents,ident,identPtrList,thisStat) + call create_IdentPtrlst(ident,ljListHead,identPtrList,thisStat) if ( thisStat /= 0 ) then status = -1 return endif + !! Allocate pointer lists allocate(point(nComponents),stat=alloc_stat) if (alloc_stat /=0) then @@ -268,13 +158,15 @@ contains ! 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 + write(default_error,*) "MPI is not set" status = -1 return endif - nrow = getNrow() - ncol = getNcol() + nrow = getNrow(plan_row) + ncol = getNcol(plan_col) !! Allocate temperary arrays to hold gather information allocate(identRow(nrow),stat=alloc_stat) if (alloc_stat /= 0 ) then @@ -287,18 +179,22 @@ contains status = -1 return endif + !! Gather idents into row and column idents + call gather(ident,identRow,plan_row) call gather(ident,identCol,plan_col) + !! Create row and col pointer lists - call new_ljatypePtrList(nrow,identRow,identPtrListRow,thisStat) + + call create_IdentPtrlst(identRow,ljListHead,identPtrListRow,thisStat) if (thisStat /= 0 ) then status = -1 return endif - call new_ljatypePtrList(ncol,identCol,identPtrListCol,thisStat) + call create_IdentPtrlst(identCol,ljListHead,identPtrListColumn,thisStat) if (thisStat /= 0 ) then status = -1 return @@ -308,69 +204,7 @@ contains deallocate(identCol) deallocate(identRow) -!! Allocate Simulation arrays -!! NOTE: This bit of code should be fixed, it can cause large -!! memory fragmentation if call repeatedly - if (.not.allocated(qRow)) then - allocate(qRow(3,nrow),stat=alloc_stat) - if (alloc_stat /= 0 ) then - status = -1 - return - endif - else - deallocate(qrow) - allocate(qRow(3,nrow),stat=alloc_stat) - if (alloc_stat /= 0 ) then - status = -1 - return - endif - endif - - if (.not.allocated(3,qCol)) then - allocate(qCol(ncol),stat=alloc_stat) - if (alloc_stat /= 0 ) then - status = -1 - return - endif - else - deallocate(qCol) - allocate(qCol(3,ncol),stat=alloc_stat) - if (alloc_stat /= 0 ) then - status = -1 - return - endif - endif - - if (.not.allocated(fRow)) then - allocate(fRow(3,nrow),stat=alloc_stat) - if (alloc_stat /= 0 ) then - status = -1 - return - endif - else - deallocate(fRow) - allocate(fRow(3,nrow),stat=alloc_stat) - if (alloc_stat /= 0 ) then - status = -1 - return - endif - endif - - if (.not.allocated(fCol)) then - allocate(fCol(3,ncol),stat=alloc_stat) - if (alloc_stat /= 0 ) then - status = -1 - return - endif - else - deallocate(fCol) - allocate(fCol(3,ncol),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,100 +236,148 @@ contains #endif - + call createMixingList(thisStat) + if (thisStat /= 0) then + status = -1 + return + endif isljFFinit = .true. - 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 + + + + tmpPtrOuter => ljListHead + tmpPtrInner => tmpPtrOuter%next + do while (associated(tmpPtrOuter)) + outerCounter = outerCounter + 1 +! 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)) + + 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 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 +#ifdef IS_MPI + real( kind = DP ), dimension(3,getNcol(plan_col)) :: 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(plan_row)) :: qRow + real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol + + real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow + real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol + real(kind = dp), dimension(3,getNlocal()) :: fMPITemp + + real(kind = dp), dimension(getNrow(plan_row)) :: eRow + real(kind = dp), dimension(getNcol(plan_col)) :: 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 @@ -504,27 +386,55 @@ 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) + if (firstTime) then + update_nlist = .true. + firstTime = .false. + endif !--------------WARNING........................... ! Zero variables, NOTE:::: Forces are zeroed in C @@ -532,44 +442,41 @@ contains ! 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 + fRow = 0.0E0_DP + fCol = 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,qRow,plan_row3) - call gather(q,qCol,plan_col3) +#ifdef IS_MPI + call gather(q,qRow,plan_row3d) + call gather(q,qCol,plan_col3d) #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 - + do i = 1, nrow point(i) = nlist + 1 -#ifdef MPI +#ifdef IS_MPI ljAtype_i => identPtrListRow(i)%this tag_i = tagRow(i) rxi = qRow(1,i) @@ -584,10 +491,10 @@ contains #endif inner: do j = j_start, ncol -#ifdef MPI +#ifdef IS_MPI ! Assign identity pointers and tags ljAtype_j => identPtrListColumn(j)%this - tag_j = tagCol(j) + tag_j = tagColumn(j) if (newtons_thrd) then if (tag_i <= tag_j) then if (mod(tag_i + tag_j,2) == 0) cycle inner @@ -608,32 +515,33 @@ contains #endif rijsq = rxij*rxij + ryij*ryij + rzij*rzij -#ifdef MPI - if (rijsq <= rlstsq .AND. & +#ifdef IS_MPI + 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 - + if (rijsq < rcutsq) then r = dsqrt(rijsq) 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 +#ifdef IS_MPI + eRow(i) = eRow(i) + pot*0.5 + eCol(i) = eCol(i) + pot*0.5 #else - pe = pe + pot + pe = pe + pot #endif efr(1,j) = -rxij @@ -647,7 +555,7 @@ contains ftmp = dudr * drdx1 -#ifdef MPI +#ifdef IS_MPI fCol(dim,j) = fCol(dim,j) - ftmp fRow(dim,i) = fRow(dim,i) + ftmp #else @@ -662,7 +570,7 @@ contains enddo inner enddo -#ifdef MPI +#ifdef IS_MPI point(nrow + 1) = nlist + 1 #else point(natoms) = nlist + 1 @@ -676,7 +584,7 @@ contains JEND = POINT(i+1) - 1 ! check thiat molecule i has neighbors if (jbeg .le. jend) then -#ifdef MPI +#ifdef IS_MPI ljAtype_i => identPtrListRow(i)%this rxi = qRow(1,i) ryi = qRow(2,i) @@ -689,11 +597,11 @@ contains #endif do jnab = jbeg, jend j = list(jnab) -#ifdef MPI +#ifdef IS_MPI 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) + 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) @@ -707,11 +615,11 @@ contains r = dsqrt(rijsq) 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 +#ifdef IS_MPI + 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 @@ -723,7 +631,7 @@ contains drdx1 = efr(dim,j) / r ftmp = dudr * drdx1 -#ifdef MPI +#ifdef IS_MPI fCol(dim,j) = fCol(dim,j) - ftmp fRow(dim,i) = fRow(dim,i) + ftmp #else @@ -739,39 +647,41 @@ contains -#ifdef MPI +#ifdef IS_MPI !!distribute forces - call scatter(fRow,f,plan_row3) + call scatter(fRow,f,plan_row3d) - call scatter(fCol,f_tmp,plan_col3) + call scatter(fCol,fMPITemp,plan_col3d) + 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) + fMPITemp(1:3,i) end do if (do_pot) then ! scatter/gather pot_row into the members of my column - call scatter(e_row,e_tmp,plan_row) + call scatter(eRow,eTemp,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) + pe_local = pe_local + eTemp(i) enddo if (newtons_thrd) then - e_tmp = 0.0E0_DP - call scatter(e_col,e_tmp,plan_col) + eTemp = 0.0E0_DP + call scatter(eCol,eTemp,plan_col) do i = 1, nlocal - pot_local = pot_local + e_tmp(i) + pe_local = pe_local + eTemp(i) enddo endif endif + potE = pe_local #endif + potE = pe + end subroutine do_lj_ff @@ -789,8 +699,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 @@ -798,7 +708,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 @@ -823,13 +733,14 @@ 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 - + + call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat) r2 = r * r @@ -846,11 +757,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 @@ -879,7 +790,7 @@ contains case ("sigma") myMixParam = 0.5_dp * (param1 + param2) - case ("epslon") + case ("epsilon") myMixParam = sqrt(param1 * param2) case default status = -1