Local minimization, Lennard-Jones potential [algorithm] by 바죠

1차 미분과 함수 값을 이용한 구조 최적화 문제의 실제 예를 보여준다.


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


LJ cluster 구조 변형 기작 시뮬레이션:
http://incredible.egloos.com/4660981


      implicit none
      integer natoms,nn
      real*8 elj
      real*8, allocatable :: x(:),v(:)
      real*8, allocatable :: g(:,:)
      integer i
      real*8 tmpx,tmpy,tmpz,obj

      natoms=7
      allocate(x(3*natoms),v(3*natoms))
      allocate(g(natoms,natoms))

      do i=1,natoms
      x(3*(i-1)+1)=0.3*i
      x(3*(i-1)+2)=0.3*i**2+1.
      x(3*(i-1)+3)=0.d0
      enddo
      i=1
      x(3*(i-1)+1)=-1.0d0
      x(3*(i-1)+2)=0.0d0
      x(3*(i-1)+3)=0.0d0
      i=2
      x(3*(i-1)+1)=1.0d0
      x(3*(i-1)+2)=0.0d0
      x(3*(i-1)+3)=0.0d0
      i=3
      x(3*(i-1)+1)=0.0d0
      x(3*(i-1)+2)=-1.0d0
      x(3*(i-1)+3)=0.0d0
      i=4
      x(3*(i-1)+1)=0.0d0
      x(3*(i-1)+2)=1.0d0
      x(3*(i-1)+3)=0.0d0
      i=5
      x(3*(i-1)+1)=1.0d0
      x(3*(i-1)+2)=1.0d0
      x(3*(i-1)+3)=0.0d0
      i=6
      x(3*(i-1)+1)=-1.0d0
      x(3*(i-1)+2)=-1.0d0
      x(3*(i-1)+3)=0.0d0
      i=7
      x(3*(i-1)+1)=0.0d0
      x(3*(i-1)+2)=0.0d0
      x(3*(i-1)+3)=0.0d0
      tmpx=0.d0 ; tmpy=0.d0 ; tmpz=0.d0
      do i=1,natoms
      tmpx=tmpx+x(3*(i-1)+1)
      tmpy=tmpy+x(3*(i-1)+2)
      tmpz=tmpz+x(3*(i-1)+3)
      enddo
      tmpx=tmpx/float(natoms) ; tmpy=tmpy/float(natoms) ; tmpz=tmpz/float(natoms)
      do i=1,natoms
      x(3*(i-1)+1)=x(3*(i-1)+1)-tmpx
      x(3*(i-1)+2)=x(3*(i-1)+2)-tmpy
      x(3*(i-1)+3)=x(3*(i-1)+3)-tmpz
      enddo

      call lj(x,v,elj,natoms,g)
      write(6,*) natoms
      write(6,*) elj
      do i=1,natoms
      write(6,'(a2,2x,3f18.8)') 'Ar',x(3*(i-1)+1),x(3*(i-1)+2),x(3*(i-1)+3)
      enddo
      write(6,*)
      do i=1,natoms
      write(6,'(3f18.8)') v(3*(i-1)+1),v(3*(i-1)+2),v(3*(i-1)+3)
      enddo
      nn=3*natoms
      call lbfgs_driv(obj,x,v,nn)
      elj=obj
      tmpx=0.d0 ; tmpy=0.d0 ; tmpz=0.d0
      do i=1,natoms
      tmpx=tmpx+x(3*(i-1)+1)
      tmpy=tmpy+x(3*(i-1)+2)
      tmpz=tmpz+x(3*(i-1)+3)
      enddo
      tmpx=tmpx/float(natoms) ; tmpy=tmpy/float(natoms) ; tmpz=tmpz/float(natoms)
      do i=1,natoms
      x(3*(i-1)+1)=x(3*(i-1)+1)-tmpx
      x(3*(i-1)+2)=x(3*(i-1)+2)-tmpy
      x(3*(i-1)+3)=x(3*(i-1)+3)-tmpz
      enddo

      write(6,*) natoms
      write(6,*) elj
      do i=1,natoms
      write(6,'(a2,2x,3f18.8)') 'Ar',x(3*(i-1)+1),x(3*(i-1)+2),x(3*(i-1)+3)
      enddo
      write(6,*)
      do i=1,natoms
      write(6,'(3f18.8)') v(3*(i-1)+1),v(3*(i-1)+2),v(3*(i-1)+3)
      enddo

      deallocate(x,v)
      deallocate(g)
      stop
      end
      subroutine lj(x,v,elj,natoms,g)
      implicit none
      integer natoms
      real*8 x(3*natoms),v(3*natoms),g(natoms,natoms)
      real*8 r6, elj, dummyx, dummyy, dummyz, xmul2, dummy, dist
      integer j1, j2, j3, j4

      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

      return
      end
!
       subroutine lbfgs_driv(theta,posi,grad,nn)
!      Written by In-Ho Lee, KRISS, January 31 (2003)
       implicit none
       integer nn
       real*8 theta,posi(nn),grad(nn)
       integer n_local,m,msave
       real*8, allocatable :: diag(:),w(:)
       real*8 eps,xtol,gtol,t1,t2,stpmin,stpmax
       integer iprint(2),iflag,icall,mp,lp,j,nwork
       logical diagco
       real*8, allocatable :: gmatwork(:,:)
       integer itrelax
!
!     The driver for LBFGS must always declare LB2 as EXTERNAL
       external lb2
       common /lb3/mp,lp,gtol,stpmin,stpmax
!
      itrelax=1000
      n_local=nn
      m=5 ; msave=7 ; nwork=n_local*(2*msave+1)+2*msave
      allocate(diag(n_local),w(nwork))
      iprint(1)= 1 ; iprint(2)= 0
!
!     We do not wish to provide the diagonal matrices Hk0, and
!     therefore set DIAGCO to FALSE.
      diagco= .false. ; eps= 1.0d-4 ; xtol= 1.0d-16 ; icall=0 ; iflag=0
!
!
      allocate(gmatwork(nn/3,nn/3))
 20   continue
      call lj(posi,grad,theta,nn/3,gmatwork)
!
      write(6,'(e14.7)') theta
      call lbfgsac(n_local,m,posi,theta,grad,diagco,diag,iprint,eps,xtol,w,iflag)
      if(iflag.le.0) go to 50
      icall=icall + 1
!     We allow at most 2000 evaluations of Function and Gradient
      if(icall >  itrelax) go to 50
      go to 20
  50  continue
      deallocate(gmatwork)
      deallocate(diag,w)
       return
       end


핑백

덧글

댓글 입력 영역

최근 포토로그



MathJax