make cp src-HCD-5D
[unres.git] / source / unres / src-HCD-5D / gradient_p.F
index 75192e9..adafa53 100644 (file)
@@ -1,17 +1,24 @@
+#ifndef LBFGS
       subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
+      include 'COMMON.CONTROL'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
       include 'COMMON.FFIELD'
       include 'COMMON.MD'
+      include 'COMMON.QRESTR'
       include 'COMMON.IOUNITS'
+      integer n,nf
+      double precision ufparm
       external ufparm
       integer uiparm(1)
       double precision urparm(1)
-      dimension x(n),g(n)
+      double precision x(n),g(n)
+      integer i,j,k,ind,ind1
+      double precision f,gthetai,gphii,galphai,gomegai
 c
 c This subroutine calculates total internal coordinate gradient.
 c Depending on the number of function evaluations, either whole energy 
@@ -30,60 +37,12 @@ c     write (iout,*) 'grad 20'
       if (nf.eq.0) return
       goto 40
    30 call var_to_geom(n,x)
-      call chainbuild 
+      call chainbuild_extconf 
 c     write (iout,*) 'grad 30'
 C
-C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
-C
-   40 call cartder
-c     write (iout,*) 'grad 40'
-c     print *,'GRADIENT: nnt=',nnt,' nct=',nct,' expon=',expon
+C Transform the gradient to the gradient in angles.
 C
-C Convert the Cartesian gradient into internal-coordinate gradient.
-C
-      ind=0
-      ind1=0
-      do i=1,nres-2
-       gthetai=0.0D0
-       gphii=0.0D0
-       do j=i+1,nres-1
-          ind=ind+1
-c         ind=indmat(i,j)
-c         print *,'GRAD: i=',i,' jc=',j,' ind=',ind
-         do k=1,3
-            gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
-          enddo
-         do k=1,3
-           gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
-          enddo
-        enddo
-       do j=i+1,nres-1
-          ind1=ind1+1
-c         ind1=indmat(i,j)
-c         print *,'GRAD: i=',i,' jx=',j,' ind1=',ind1
-         do k=1,3
-           gthetai=gthetai+dxdv(k,ind1)*gradx(k,j,icg)
-           gphii=gphii+dxdv(k+3,ind1)*gradx(k,j,icg)
-          enddo
-        enddo
-       if (i.gt.1) g(i-1)=gphii
-       if (n.gt.nphi) g(nphi+i)=gthetai
-      enddo
-      if (n.le.nphi+ntheta) goto 10
-      do i=2,nres-1
-       if (itype(i).ne.10) then
-          galphai=0.0D0
-         gomegai=0.0D0
-         do k=1,3
-           galphai=galphai+dxds(k,i)*gradx(k,i,icg)
-          enddo
-         do k=1,3
-           gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
-          enddo
-          g(ialph(i,1))=galphai
-         g(ialph(i,1)+nside)=gomegai
-        endif
-      enddo
+   40 call cart2intgrad(n,g)
 C
 C Add the components corresponding to local energy terms.
 C
@@ -109,7 +68,7 @@ cd    write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n)
       end
 C-------------------------------------------------------------------------
       subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
@@ -117,10 +76,14 @@ C-------------------------------------------------------------------------
       include 'COMMON.INTERACT'
       include 'COMMON.FFIELD'
       include 'COMMON.IOUNITS'
+      integer n,nf
+      double precision ufparm
       external ufparm
       integer uiparm(1)
       double precision urparm(1)
-      dimension x(maxvar),g(maxvar)
+      double precision x(maxvar),g(maxvar),gg(maxvar)
+      integer i,j,k,ig,ind,ij,igall
+      double precision f,gthetai,gphii,galphai,gomegai
 
       icg=mod(nf,2)+1
       if (nf-nfl+1) 20,30,40
@@ -148,58 +111,33 @@ c      write(iout,*) (var(i),i=1,nvar)
 C
 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 C
-   40 call cartder
+   40 call cart2intgrad(n,gg)
 C
 C Convert the Cartesian gradient into internal-coordinate gradient.
 C
 
       ig=0
-      ind=nres-2                                                                    
+      ind=nres-2 
       do i=2,nres-2                
-       IF (mask_phi(i+2).eq.1) THEN                                             
-        gphii=0.0D0                                                             
-        do j=i+1,nres-1                                                         
-          ind=ind+1                                 
-          do k=1,3                                                              
-            gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)                            
-            gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)                           
-          enddo                                                                 
-        enddo                                                                   
+       IF (mask_phi(i+2).eq.1) THEN
         ig=ig+1
-        g(ig)=gphii
-       ELSE
-        ind=ind+nres-1-i
+        g(ig)=gg(i-1)
        ENDIF
       enddo                                        
 
 
-      ind=0
       do i=1,nres-2
        IF (mask_theta(i+2).eq.1) THEN
         ig=ig+1
-       gthetai=0.0D0
-       do j=i+1,nres-1
-          ind=ind+1
-         do k=1,3
-            gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
-            gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
-          enddo
-        enddo
-        g(ig)=gthetai
-       ELSE
-        ind=ind+nres-1-i
+        g(ig)=gg(nphi+i)
        ENDIF
       enddo
 
       do i=2,nres-1
-       if (itype(i).ne.10) then
+        if (itype(i).ne.10) then
          IF (mask_side(i).eq.1) THEN
           ig=ig+1
-          galphai=0.0D0
-         do k=1,3
-           galphai=galphai+dxds(k,i)*gradx(k,i,icg)
-          enddo
-          g(ig)=galphai
+          g(ig)=gg(ialph(i,1))
          ENDIF
         endif
       enddo
@@ -209,11 +147,7 @@ C
         if (itype(i).ne.10) then
          IF (mask_side(i).eq.1) THEN
           ig=ig+1
-         gomegai=0.0D0
-         do k=1,3
-           gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
-          enddo
-         g(ig)=gomegai
+          g(ig)=gg(ialph(i,1)+nside)
          ENDIF
         endif
       enddo
@@ -257,21 +191,25 @@ cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
 cd      enddo
       return
       end
+#endif
 C-------------------------------------------------------------------------
       subroutine cartgrad
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
 #endif
+      include 'COMMON.CONTROL'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
       include 'COMMON.FFIELD'
       include 'COMMON.MD'
+      include 'COMMON.QRESTR'
       include 'COMMON.IOUNITS'
       include 'COMMON.TIME1'
+      integer i,j,kk
 c
 c This subrouting calculates total Cartesian coordinate gradient. 
 c The subroutine chainbuild_cart and energy MUST be called beforehand.
@@ -376,9 +314,73 @@ cd      write(iout,*) 'calling int_to_cart'
 #endif
       return
       end
+c---------------------------------------------------------------------------
+#ifdef FIVEDIAG
+      subroutine grad_transform
+      implicit none
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.CONTROL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      include 'COMMON.MD'
+      include 'COMMON.QRESTR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.TIME1'
+      integer i,j,kk
+#ifdef DEBUG
+      write (iout,*)"Converting virtual-bond gradient to CA/SC gradient"
+      write (iout,*) "dC/dX gradient"
+      do i=0,nres
+        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
+     &      (gxcart(j,i),j=1,3)
+      enddo
+#endif
+      do i=nres,1,-1
+        do j=1,3
+          gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+!          gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+        enddo
+!        write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
+!            (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
+      enddo
+! Correction: dummy residues
+      do i=2,nres
+        if (itype(i-1).eq.ntyp1 .and. itype(i).ne.ntyp1) then
+          gcart(:,i)=gcart(:,i)+gcart(:,i-1)
+        else if (itype(i-1).ne.ntyp1 .and. itype(i).eq.ntyp1) then
+          gcart(:,i-1)=gcart(:,i-1)+gcart(:,i)
+        endif
+      enddo
+c      if (nnt.gt.1) then
+c        do j=1,3
+c          gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
+c        enddo
+c      endif
+c      if (nct.lt.nres) then
+c        do j=1,3
+c!          gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
+c          gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
+c        enddo
+c      endif
+#ifdef DEBUG
+      write (iout,*) "CA/SC gradient"
+      do i=1,nres
+        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
+     &      (gxcart(j,i),j=1,3)
+      enddo
+#endif
+      return
+      end
+#endif
 C-------------------------------------------------------------------------
       subroutine zerograd
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.DERIV'
       include 'COMMON.CHAIN'
@@ -386,6 +388,7 @@ C-------------------------------------------------------------------------
       include 'COMMON.MD'
       include 'COMMON.SCCOR'
       include 'COMMON.SHIELD'
+      integer i,j,kk,intertyp,maxshieldlist
       maxshieldlist=0
 C
 C Initialize Cartesian-coordinate gradient
@@ -461,14 +464,18 @@ C           grad_shield_side_ca(j,kk,i)=0.0d0
           do intertyp=1,3
            gloc_sc(intertyp,i,icg)=0.0d0
           enddo
+        enddo
+      enddo
 #ifndef DFA
+      do i=1,nres
+        do j=1,3
           gdfad(j,i)=0.0d0
           gdfat(j,i)=0.0d0
           gdfan(j,i)=0.0d0
           gdfab(j,i)=0.0d0
-#endif
         enddo
       enddo
+#endif
 C
 C Initialize the gradient of local energy terms.
 C