string method (LJ potential) [algorithm] by 바죠

string method (LJ potential) [algorithm]
cf. Nudged Elastic Band (NEB) method

transition pathway/ minimum energy path/ transition state를 찾는 방법 중 하나로 알려진 string method의 간단한 사용예를 보였다.

참고:
http://web.math.princeton.edu/string/

http://www.math.nus.edu.sg/~matrw/string/codes/ZTS.html

LJ potential:
http://incredible.egloos.com/4242358

LJ cluster 구조 최적화:
http://incredible.egloos.com/4659953

!234567890
       module string
!      Written by In-Ho Lee, KRISS, December 22 (2011); based on  W. Ren's string programs in fortran 77.
       implicit none
       private
       save
       integer md,n
       real*8, allocatable :: xxx(:,:),eee(:)
       real*8, allocatable :: r(:),t(:,:),s(:),s1(:),x1(:,:),e(:)
       real*8, allocatable :: gwork(:,:),fwork(:)
       public :: md,n,xxx,eee
       public :: initial,final,mvtan,mvper,print_intermediate

       contains

       subroutine initial(dt_dum)
       implicit none
       real*8 dt_dum
       integer i,j,natoms
       real*8 u,tmpx,tmpy,tmpz
       logical lexist
       character*2 ichar

       dt_dum=1.d-5
       dt_dum=1.d-6
       dt_dum=1.d-7
       md=3*7
       n=30
       n=20
       write(6,*) md,n,' md,n,dt ',dt_dum
       allocate(xxx(md,0:n))
       allocate(eee(0:n))
       allocate(r(md),t(md,n-1))
       allocate(s(0:n),s1(0:n),x1(md,0:n))
       allocate(e(0:n))
       allocate(gwork(md/3,md/3))
       allocate(fwork(md))

       inquire(file='confi.a',exist=lexist)
       if(.not. lexist)then
       write(6,*) 'file confi.a is absent'
                       stop
                       endif
       open(1,file='confi.a',form='formatted')
       read(1,*) natoms
       if(md/3 /= natoms)then
       write(6,*) 'check confi.a'
                         stop
                         endif
       read(1,*)
       do i=1,natoms
       read(1,*) ichar,xxx(3*(i-1)+1,0),xxx(3*(i-1)+2,0),xxx(3*(i-1)+3,0)
       enddo
       close(1)
       j=0
       call shift_conformations(j)

       inquire(file='confi.b',exist=lexist)
       if(.not. lexist)then
       write(6,*) 'file confi.b is absent'
                       stop
                       endif
       open(2,file='confi.b',form='formatted')
       read(2,*) i
       if(md/3 /= i     )then
       write(6,*) 'check confi.b'
                         stop
                         endif
       read(2,*)
       do i=1,natoms
       read(2,*) ichar, xxx(3*(i-1)+1,n),xxx(3*(i-1)+2,n),xxx(3*(i-1)+3,n)
       enddo
       close(2)
       j=n
       call shift_conformations(j)

       do j=1,n-1
       u=dble(j)/dble(n)
       do i=1,md
       xxx(i,j)=xxx(i,0)*(1.d0-u)+xxx(i,n)*u
       enddo
       enddo

       inquire(file='fort.8',exist=lexist)
       if(lexist)then
       write(6,*) 'file fort.8 is present'
       write(6,*) 'iterative calculation'
       open(8,file='fort.8',form='formatted')
       do j=0,n
       read(8,*) i
       if(md/3 /= i)then
       write(6,*) 'check fort.8'
                    stop
                    endif
       read(8,*)
       do i=1,md/3
       read(8,*) ichar,xxx(3*(i-1)+1,j),xxx(3*(i-1)+2,j),xxx(3*(i-1)+3,j)
       enddo
       enddo
       close(8)
       do j=0,n
       call shift_conformations(j)
       enddo
                 endif
       end subroutine initial

       subroutine final
       implicit none

       call print_intermediate
       write(6,*) 'cp fort.8 test.xyz ; xmakemol -f test.xyz'
       write(6,*) 'xmgrace fort.34'

       deallocate(xxx)
       deallocate(eee)
       deallocate(r,t)
       deallocate(s,s1,x1)
       deallocate(e)
       deallocate(gwork)
       deallocate(fwork)
       end subroutine final

       subroutine print_intermediate
       implicit none
       integer j,i

       open(34,file='fort.34',form='formatted')
       open(8,file='fort.8',form='formatted')
       do j=0,n
       call shift_conformations(j)
       write(8,*) md/3
       write(8,*) eee(j)
       write(34,*) j,eee(j)
       do i=1,md/3
       write(8,'(a2,3f18.8)') 'Ar', xxx(3*(i-1)+1,j),xxx(3*(i-1)+2,j),xxx(3*(i-1)+3,j)
       enddo
       enddo
       close(8)
       close(34)
       end subroutine print_intermediate

       subroutine shift_conformations(j)
       implicit none
       integer j
       real*8 tmpx,tmpy,tmpz
       integer i,natoms

       natoms=md/3
       tmpx=0.d0 ; tmpy=0.d0 ; tmpz=0.d0
       do i=1,md/3
       tmpx=tmpx+xxx(3*(i-1)+1,j)
       tmpy=tmpy+xxx(3*(i-1)+2,j)
       tmpz=tmpz+xxx(3*(i-1)+3,j)
       enddo
       tmpx=tmpx/float(natoms) ; tmpy=tmpy/float(natoms) ; tmpz=tmpz/float(natoms)
       do i=1,natoms
       xxx(3*(i-1)+1,j)=xxx(3*(i-1)+1,j)-tmpx
       xxx(3*(i-1)+2,j)=xxx(3*(i-1)+2,j)-tmpy
       xxx(3*(i-1)+3,j)=xxx(3*(i-1)+3,j)-tmpz
       enddo
       end subroutine shift_conformations

      subroutine mvper(dt,x,error)
      implicit none
! Evlove each image by projected potential force
! by W. Ren
!
! input     dt: time step
!            x: coordinates of images along the string
! Output     x: new position of the string after one time step
!        error: the L_\inf norm of the projected force
!
      real*8 dt,error,x(md,0:n)
      integer k,i
      real*8 err,err1,tmp,s
!!    real*8 r(md),t(md,n-1)
!
! calculate the unit tangent vector along the string
      call tangent(x,t)
!
      error=0.d0
      do 1 k=1,n-1
!
! calculate the potential force of the k-th image
        call force(md,x(1,k),r)
!
! project the force onto the plane normal to the string
        s=0.d0
        do 4 i=1,md
          s=s+r(i)*t(i,k)
  4     continue
        do 2 i=1,md
          r(i)=r(i)-s*t(i,k)
  2     continue
!
! Update the images
        tmp=0.d0
        do 3 i=1,md
          x(i,k)=x(i,k)+dt*r(i)
          tmp=tmp+r(i)**2.d0
  3     continue
        tmp=dsqrt(tmp)
        error=dmax1(error,tmp)
!
  1   continue
      end subroutine mvper

      subroutine mvtan(x)
      implicit none
! Reparameterization of the string by polynomial interpolation
! By W. Ren
!
      real*8 x(md,0:n)
      integer k,i
!!    real*8 s(0:n),s1(0:n),x1(md,0:n)
!
! calculate the distance between neighboring images along the string
      do 1 k=1,n
        s(k)=0.d0
        do 2 i=1,md
          s(k)=s(k)+(x(i,k)-x(i,k-1))**2.d0
  2     continue
        s(k)=dsqrt(s(k))
  1   continue
!
! normaized arclength starting from the initial state
      s(0)=0.d0
      do 3 k=1,n
        s(k)=s(k-1)+s(k)
  3   continue
      do 31 k=1,n
        s(k)=s(k)/s(n)
  31  continue
!
! Equal-arclength: the locations at which new images are calculated
! by interpolation
      do 4 k=0,n
        s1(k)=dble(k)/dble(n)
  4   continue
!
! interpolation to get new images
!
      call intpol(md,n,s,x,n,s1,x1)
!
      call dcopy(md*(n+1),x1,1,x,1)
      end subroutine mvtan

      subroutine tangent(x,t)
      implicit none
! calculate the unit tangent vectors along the string
! by W. Ren
!
! Input  x: coordinates of the images along the string
! Output t: unit tangent vectors
!
      real*8 x(md,0:n),t(md,n-1)
      real*8 e1,e2,tmp1,tmp2
      integer k,i
!!    real*8 e(0:n)
!
! upwind scheme based on the potential energy
!
      call energy(x,e)
!
      do 1 k=1,n-1
        if(e(k+1).ge.e(k).and.e(k).ge.e(k-1)) then
          do 2 i=1,md
            t(i,k)=x(i,k+1)-x(i,k)
  2       continue
        else if(e(k+1).le.e(k).and.e(k).le.e(k-1)) then
          do 3 i=1,md
            t(i,k)=x(i,k)-x(i,k-1)
  3       continue
        else
          tmp1=dabs(e(k-1)-e(k))
          tmp2=dabs(e(k)-e(k+1))
          e1=dmax1(tmp1,tmp2)
          e2=dmin1(tmp1,tmp2)
          if(e(k+1).ge.e(k-1)) then
            do 4 i=1,md
              t(i,k)=e1*(x(i,k+1)-x(i,k))+    e2*(x(i,k)-x(i,k-1))
  4         continue
          else
            do 5 i=1,md
              t(i,k)=e2*(x(i,k+1)-x(i,k))    +e1*(x(i,k)-x(i,k-1))
  5         continue
          endif
        endif
! nomalization
        tmp1=0.d0
        do 6 i=1,md
          tmp1=tmp1+t(i,k)**2.d0
  6     continue
        tmp1=dsqrt(tmp1)
        do 7 i=1,md
          t(i,k)=t(i,k)/tmp1
  7     continue
!
  1   continue
      end subroutine tangent

      subroutine energy(xyz,e)
      implicit none
! Calculate the potential (Mueller) energy of each image along the path
! by W. Ren
!
! Input  : xyz(2,0:n), the coordinates of each image along the string
! Output : e(0:n), the potential energy of each image
!
      real*8 xyz(md,0:n),e(0:n)
      integer k
      real*8 t1,t2,t3,t4,t5,t6,f,g,f1,g1
!
      do k=0,n
      call lj(xyz(1,k),fwork,e(k),md/3,gwork)
      enddo
      eee=e
      end subroutine energy

      subroutine force(md,xyz,ff)
      implicit none
      integer md
      real*8 xyz(md),ff(md)
      real*8 tmp
      integer i
! Calculate the potential (Mueller) force
! by W. Ren
!
! Input xyz: coordinate of the image
! output ff: potential force
!
      call lj(xyz,ff,tmp,md/3,gwork)
      do i=1,md
      ff(i)=-ff(i)
      enddo
      end subroutine force

       end module string

      subroutine lj(x,v,elj,natoms,g)
!      Written by In-Ho Lee, KRISS, December 22 (2011); based on  Wales's GMIN programs in fortran 77.
      implicit none
      integer natoms
      real*8 x(3*natoms),v(3*natoms),elj,g(natoms,natoms)
      integer j1, j2, j3, j4
      real*8 r6, dummyx, dummyy, dummyz, xmul2, dummy, dist

      elj=0.0d0
         do j1=1,natoms
            j3=3*j1
            g(j1,j1)=0.0d0
            do j2=j1+1,natoms
               j4=3*j2
               dist=(x(j3-2)-x(j4-2))**2+(x(j3-1)-x(j4-1))**2+(x(j3)-x(j4))**2
               dist=1.0d0/dist
               r6=dist**3
               dummy=r6*(r6-1.0d0)
               elj=elj+dummy
               dist=dist*r6
               g(j2,j1)=-24.0d0*(2.0d0*r6-1.0d0)*dist
               g(j1,j2)=g(j2,j1)
            enddo
         enddo
      elj=elj*4.0d0

         do j1=1,natoms
            j3=3*j1
            dummyx=0.0d0
            dummyy=0.0d0
            dummyz=0.0d0
            do j4=1,natoms
               j2=3*j4
               xmul2=g(j4,j1)
               dummyx=dummyx+xmul2*(x(j3-2)-x(j2-2))
               dummyy=dummyy+xmul2*(x(j3-1)-x(j2-1))
               dummyz=dummyz+xmul2*(x(j3)  -x(j2))
            enddo
            v(j3-2)=dummyx
            v(j3-1)=dummyy
            v(j3)=dummyz
         enddo
      end



       program string_driver
!      Written by In-Ho Lee, KRISS, December 22 (2011); based on  W. Ren's string programs in fortran 77.
       USE string, ONLY : md,n,xxx,eee
       USE string, ONLY : initial,final,mvtan,mvper,print_intermediate
       implicit none
       real*8 dt,error
       integer ic
       character*2 ichar

       dt=0.0004d0
       dt=1.0d-7
       call initial(dt)
!
       ic=0
  100  continue
       ic=ic+1
       if(mod(ic,10).eq.1) then
       call mvtan(xxx)
                           endif
       call mvper(dt,xxx,error)

       if(mod(ic,10).eq.1) then
       write(*,*) 'ic=',ic, '  error =',error
                           endif

       if(mod(ic,1000).eq.0) then
       call print_intermediate
                             endif

       if(error.lt.1.e-6) then
                          goto 101
                          else
       if(ic > 1000000)then
       write(6,*) 'ic=',ic, '  error =',error
                       goto 101
                       endif
                          endif
       goto 100
  101  continue

       call final
!
       stop
       end program string_driver





핑백

덧글

댓글 입력 영역

최근 포토로그



MathJax