--- trunk/mdtools/md_code/lj_FF.F90 2003/01/24 21:36:52 246 +++ trunk/mdtools/md_code/lj_FF.F90 2003/01/27 18:28:11 247 @@ -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.10 2003-01-24 21:36:52 chuckv Exp $, $Date: 2003-01-24 21:36:52 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $ +!! @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 + 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 ) :: epsilon = 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(:,:), pointer :: ljMixed => null() -!! identity pointer list for force loop. - type (lj_atypePtr), dimension(:), pointer :: identPtrList => null() !! Neighbor list and commom arrays @@ -64,29 +36,15 @@ 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. -!! 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 - -!! Row potential energy array - real( kind = dp ), allocatable, dimension(:) :: eRow -!! Column potential energy array - real( kind = dp ), allocatable, dimension(:) :: eColumn - - +!! 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 @@ -105,6 +63,7 @@ contains contains +!! 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) :: epsilon @@ -112,143 +71,41 @@ contains 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), dimension(:,:), pointer :: thisMix - type (lj_atype), 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%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 + 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%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(ljMixed) - 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 - - 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 - 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)) 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 - ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma = & - calcLJMix("sigma",this_lj_atype%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 - - 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 - - ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%epsilon = & - calcLJMix("epsilon",this_lj_atype%epsilon, & - lj_atype_ptr%epsilon) - -! Advance to next pointer - 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 - -! Increment number of atypes - n_lj_atypes = n_lj_atypes + 1 -! Set up self mixing end subroutine new_lj_atype @@ -274,16 +131,14 @@ contains integer :: alloc_stat #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) + call create_IdentPtrlst(ident,ljListHead,identPtrList,thisStat) if ( thisStat /= 0 ) then status = -1 return endif - write(*,*) "component pointer list initialized" + !! Allocate pointer lists allocate(point(nComponents),stat=alloc_stat) if (alloc_stat /=0) then @@ -323,13 +178,13 @@ contains call gather(ident,identCol,plan_col) !! Create row and col pointer lists - call new_ljatypePtrList(nrow,identRow,identPtrListRow,thisStat) + call create_IdentPtrlst(identRow,ljListTail,identPtrListRow,thisStat) if (thisStat /= 0 ) then status = -1 return endif - call new_ljatypePtrList(ncol,identCol,identPtrListCol,thisStat) + call create_IdentPtrlst(identCol,ljListTail,identPtrListCol,thisStat) if (thisStat /= 0 ) then status = -1 return @@ -339,116 +194,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 - - 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) @@ -480,11 +226,11 @@ contains #endif - - do i = 1, nComponents - write(*,*) "Ident pointer list ident: ", identPtrList(i)%this%atype_ident - end do - + call createMixingList(thisStat) + if (thisStat /= 0) then + status = -1 + return + endif isljFFinit = .true. end subroutine init_ljFF @@ -493,82 +239,99 @@ contains -!! 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,dimension(:), intent(in) :: ident - integer, optional :: status - type(lj_atypePtr), dimension(:), pointer :: 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 - write(*,*) "Creating new ljatypePtrList...." -! First time through, allocate list - 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(PtrList) - allocate(PtrList(mysize)) - endif + listSize = getListLen(ljListHead) + if (listSize == 0) then + status = -1 + return + end if + -! 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 + 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), 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( 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) @@ -589,6 +352,22 @@ contains 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