make cp src-HCD-5D
[unres.git] / source / unres / src-HCD-5D / checkder_p.F
index 03df287..48eedda 100644 (file)
@@ -1,181 +1,7 @@
-      subroutine check_cartgrad
-C Check the gradient of Cartesian coordinates in internal coordinates.
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'COMMON.CONTROL'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.VAR'
-      include 'COMMON.CHAIN'
-      include 'COMMON.GEO'
-      include 'COMMON.LOCAL'
-      include 'COMMON.DERIV'
-      dimension temp(6,maxres),xx(3),gg(3)
-      indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
-*
-* Check the gradient of the virtual-bond and SC vectors in the internal
-* coordinates.
-*    
-      print '("Calling CHECK_ECART",1pd12.3)',aincr
-      write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr
-      aincr2=0.5d0*aincr
-      call cartder
-      write (iout,'(a)') '**************** dx/dalpha'
-      write (iout,'(a)')
-      do i=2,nres-1
-       alphi=alph(i)
-       alph(i)=alph(i)+aincr
-       do k=1,3
-         temp(k,i)=dc(k,nres+i)
-        enddo
-       call chainbuild
-       do k=1,3
-         gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
-         xx(k)=dabs((gg(k)-dxds(k,i))/(aincr*dabs(dxds(k,i))+aincr))
-        enddo
-        write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)') 
-     &  i,(gg(k),k=1,3),(dxds(k,i),k=1,3),(xx(k),k=1,3)
-        write (iout,'(a)')
-       alph(i)=alphi
-       call chainbuild
-      enddo
-      write (iout,'(a)')
-      write (iout,'(a)') '**************** dx/domega'
-      write (iout,'(a)')
-      do i=2,nres-1
-       omegi=omeg(i)
-       omeg(i)=omeg(i)+aincr
-       do k=1,3
-         temp(k,i)=dc(k,nres+i)
-        enddo
-       call chainbuild
-       do k=1,3
-          gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
-          xx(k)=dabs((gg(k)-dxds(k+3,i))/
-     &          (aincr*dabs(dxds(k+3,i))+aincr))
-        enddo
-        write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)') 
-     &      i,(gg(k),k=1,3),(dxds(k+3,i),k=1,3),(xx(k),k=1,3)
-        write (iout,'(a)')
-       omeg(i)=omegi
-       call chainbuild
-      enddo
-      write (iout,'(a)')
-      write (iout,'(a)') '**************** dx/dtheta'
-      write (iout,'(a)')
-      do i=3,nres
-       theti=theta(i)
-        theta(i)=theta(i)+aincr
-        do j=i-1,nres-1
-          do k=1,3
-            temp(k,j)=dc(k,nres+j)
-          enddo
-        enddo
-        call chainbuild
-        do j=i-1,nres-1
-         ii = indmat(i-2,j)
-c         print *,'i=',i-2,' j=',j-1,' ii=',ii
-         do k=1,3
-           gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
-           xx(k)=dabs((gg(k)-dxdv(k,ii))/
-     &            (aincr*dabs(dxdv(k,ii))+aincr))
-          enddo
-          write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') 
-     &        i,j,(gg(k),k=1,3),(dxdv(k,ii),k=1,3),(xx(k),k=1,3)
-          write(iout,'(a)')
-        enddo
-        write (iout,'(a)')
-        theta(i)=theti
-        call chainbuild
-      enddo
-      write (iout,'(a)') '***************** dx/dphi'
-      write (iout,'(a)')
-      do i=4,nres
-        phi(i)=phi(i)+aincr
-        do j=i-1,nres-1
-          do k=1,3
-            temp(k,j)=dc(k,nres+j)
-          enddo
-        enddo
-        call chainbuild
-        do j=i-1,nres-1
-         ii = indmat(i-2,j)
-c         print *,'ii=',ii
-         do k=1,3
-           gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
-            xx(k)=dabs((gg(k)-dxdv(k+3,ii))/
-     &            (aincr*dabs(dxdv(k+3,ii))+aincr))
-          enddo
-          write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') 
-     &        i,j,(gg(k),k=1,3),(dxdv(k+3,ii),k=1,3),(xx(k),k=1,3)
-          write(iout,'(a)')
-        enddo
-        phi(i)=phi(i)-aincr
-        call chainbuild
-      enddo
-      write (iout,'(a)') '****************** ddc/dtheta'
-      do i=1,nres-2
-        thet=theta(i+2)
-        theta(i+2)=thet+aincr
-        do j=i,nres
-          do k=1,3 
-            temp(k,j)=dc(k,j)
-          enddo
-        enddo
-        call chainbuild 
-        do j=i+1,nres-1
-         ii = indmat(i,j)
-c         print *,'ii=',ii
-         do k=1,3
-           gg(k)=(dc(k,j)-temp(k,j))/aincr
-           xx(k)=dabs((gg(k)-dcdv(k,ii))/
-     &           (aincr*dabs(dcdv(k,ii))+aincr))
-          enddo
-          write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') 
-     &           i,j,(gg(k),k=1,3),(dcdv(k,ii),k=1,3),(xx(k),k=1,3)
-         write (iout,'(a)')
-        enddo
-        do j=1,nres
-          do k=1,3
-            dc(k,j)=temp(k,j)
-          enddo 
-        enddo
-        theta(i+2)=thet
-      enddo    
-      write (iout,'(a)') '******************* ddc/dphi'
-      do i=1,nres-3
-        phii=phi(i+3)
-        phi(i+3)=phii+aincr
-        do j=1,nres
-          do k=1,3 
-            temp(k,j)=dc(k,j)
-          enddo
-        enddo
-        call chainbuild 
-        do j=i+2,nres-1
-         ii = indmat(i+1,j)
-c         print *,'ii=',ii
-         do k=1,3
-           gg(k)=(dc(k,j)-temp(k,j))/aincr
-            xx(k)=dabs((gg(k)-dcdv(k+3,ii))/
-     &           (aincr*dabs(dcdv(k+3,ii))+aincr))
-          enddo
-          write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') 
-     &         i,j,(gg(k),k=1,3),(dcdv(k+3,ii),k=1,3),(xx(k),k=1,3)
-         write (iout,'(a)')
-        enddo
-        do j=1,nres
-          do k=1,3
-            dc(k,j)=temp(k,j)
-          enddo
-        enddo
-        phi(i+3)=phii   
-      enddo   
-      return
-      end
 C----------------------------------------------------------------------------
       subroutine check_ecart
 C Check the gradient of the energy in Cartesian coordinates. 
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.CONTROL'
       include 'COMMON.CHAIN'
@@ -183,12 +9,20 @@ C Check the gradient of the energy in Cartesian coordinates.
       include 'COMMON.IOUNITS'
       include 'COMMON.VAR'
       include 'COMMON.CONTACTS'
+      integer i,j,k
+      integer icall
       common /srutu/ icall
-      dimension ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),g(maxvar)
-      dimension grad_s(6,maxres)
+      double precision ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
+     & g(maxvar),grad_s(6,maxres)
       double precision energia(0:n_ene),energia1(0:n_ene)
+      double precision aincr2,etot,etot1,etot2
+      double precision dist,alpha,beta
+      double precision funcgrad,ff
+      external funcgrad
+      integer nf
       integer uiparm(1)
       double precision urparm(1)
+      double precision fdum
       external fdum
       icg=1
       nf=0
@@ -202,7 +36,11 @@ C Check the gradient of the energy in Cartesian coordinates.
       call etotal(energia(0))
       etot=energia(0)
       call enerprint(energia(0))
+#ifdef LBFGS
+      ff=funcgrad(x,g)
+#else
       call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
+#endif
       icall =1
       do i=1,nres
         write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
@@ -253,7 +91,7 @@ C Check the gradient of the energy in Cartesian coordinates.
 c----------------------------------------------------------------------------
       subroutine check_ecartint
 C Check the gradient of the energy in Cartesian coordinates. 
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.CONTROL'
       include 'COMMON.CHAIN'
@@ -264,23 +102,27 @@ C Check the gradient of the energy in Cartesian coordinates.
       include 'COMMON.MD'
       include 'COMMON.LOCAL'
       include 'COMMON.SPLITELE'
+      integer icall
       common /srutu/ icall
-      dimension ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
-     &  g(maxvar)
-      dimension dcnorm_safe(3),dxnorm_safe(3)
-      dimension grad_s(6,0:maxres),grad_s1(6,0:maxres)
+      double precision ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),
+     & x(maxvar),g(maxvar)
+      double precision dcnorm_safe(3),dxnorm_safe(3)
+      double precision grad_s(6,0:maxres),grad_s1(6,0:maxres)
       double precision phi_temp(maxres),theta_temp(maxres),
      &  alph_temp(maxres),omeg_temp(maxres)
       double precision energia(0:n_ene),energia1(0:n_ene)
       integer uiparm(1)
       double precision urparm(1)
       external fdum
+      integer i,j,k,nf
+      double precision etot,etot1,etot2,etot11,etot12,etot21,etot22
+      double precision dist,alpha,beta
 c      r_cut=2.0d0
 c      rlambd=0.3d0
       icg=1
       nf=0
       nfl=0                
-      print *,"ATU 3"
+c      print *,"ATU 3"
       call int_from_cart1(.false.)
       call intout
 c      call intcartderiv
@@ -325,15 +167,15 @@ c        call flush(iout)
         call etotal_long(energia(0))
         call enerprint(energia(0))
         call flush(iout)
-        write (iout,*) "enter cartgrad"
-        call flush(iout)
+c        write (iout,*) "enter cartgrad"
+c        call flush(iout)
         call cartgrad
-        write (iout,*) "exit cartgrad"
-        call flush(iout)
+c        write (iout,*) "exit cartgrad"
+c        call flush(iout)
         icall =1
         write (iout,*) "longrange grad"
         do i=1,nres
-          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
+          write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
      &    (gxcart(j,i),j=1,3)
         enddo
         do j=1,3
@@ -349,15 +191,15 @@ c        call flush(iout)
         call etotal_short(energia(0))
         call enerprint(energia(0))
         call flush(iout)
-        write (iout,*) "enter cartgrad"
-        call flush(iout)
+c        write (iout,*) "enter cartgrad"
+c        call flush(iout)
         call cartgrad
-        write (iout,*) "exit cartgrad"
-        call flush(iout)
+c        write (iout,*) "exit cartgrad"
+c        call flush(iout)
         icall =1
         write (iout,*) "shortrange grad"
         do i=1,nres
-          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
+          write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
      &    (gxcart(j,i),j=1,3)
         enddo
         do j=1,3
@@ -493,7 +335,7 @@ c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
       end
 c-------------------------------------------------------------------------
       subroutine int_from_cart1(lprn)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
@@ -509,6 +351,10 @@ c-------------------------------------------------------------------------
       include 'COMMON.SETUP'
       include 'COMMON.TIME1'
       logical lprn 
+      integer i,j
+      double precision dnorm1,dnorm2,be
+      double precision time00
+      double precision dist,alpha,beta
       if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
 #ifdef TIMING
       time01=MPI_Wtime()
@@ -635,7 +481,7 @@ cd       call flush(iout)
 c----------------------------------------------------------------------------
       subroutine check_eint
 C Check the gradient of energy in internal coordinates.
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.CONTROL'
       include 'COMMON.CHAIN'
@@ -643,14 +489,20 @@ C Check the gradient of energy in internal coordinates.
       include 'COMMON.IOUNITS'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
+      integer icall
       common /srutu/ icall
-      dimension x(maxvar),gana(maxvar),gg(maxvar)
+      double precision x(maxvar),gana(maxvar),gg(maxvar)
       integer uiparm(1)
       double precision urparm(1)
       double precision energia(0:n_ene),energia1(0:n_ene),
      &  energia2(0:n_ene)
       character*6 key
+      double precision fdum
       external fdum
+      double precision funcgrad,ff
+      external funcgrad
+      integer i,ii,nf
+      double precision xi,etot,etot1,etot2
       call zerograd
 c      aincr=1.0D-7
       print '("Calling CHECK_INT",1pd12.3)',aincr
@@ -678,7 +530,15 @@ c      aincr=1.0D-7
       nf=1
       nfl=3
 cd    write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
+c      write (iout,*) "Before gradient"
+c      call flush(iout)
+#ifdef LBFGS
+      ff=funcgrad(x,gana)
+#else
       call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
+#endif
+c      write (iout,*) "After gradient"
+c      call flush(iout)
 cd    write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
       icall=1
       do i=1,nvar
@@ -694,7 +554,7 @@ cd    write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
         call etotal(energia2(0))
         etot2=energia2(0)
         gg(i)=(etot2-etot1)/aincr
-        write (iout,*) i,etot1,etot2
+c        write (iout,*) i,etot1,etot2
         x(i)=xi
       enddo
       write (iout,'(/2a)')' Variable        Numerical       Analytical',