Merge branch 'lipid' of mmka.chem.univ.gda.pl:unres into lipid
[unres.git] / source / unres / src_MD-NEWSC / gradient_p.F
diff --git a/source/unres/src_MD-NEWSC/gradient_p.F b/source/unres/src_MD-NEWSC/gradient_p.F
new file mode 100644 (file)
index 0000000..7fec1e8
--- /dev/null
@@ -0,0 +1,421 @@
+      subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      include 'COMMON.MD'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.SCCOR'
+      external ufparm
+      integer uiparm(1)
+      double precision urparm(1)
+      dimension x(maxvar),g(maxvar)
+c
+c This subroutine calculates total internal coordinate gradient.
+c Depending on the number of function evaluations, either whole energy 
+c is evaluated beforehand, Cartesian coordinates and their derivatives in 
+c internal coordinates are reevaluated or only the cartesian-in-internal
+c coordinate derivatives are evaluated. The subroutine was designed to work
+c with SUMSL.
+c 
+c
+      icg=mod(nf,2)+1
+
+cd      print *,'grad',nf,icg
+      if (nf-nfl+1) 20,30,40
+   20 call func(n,x,nf,f,uiparm,urparm,ufparm)
+c     write (iout,*) 'grad 20'
+      if (nf.eq.0) return
+      goto 40
+   30 call var_to_geom(n,x)
+      call chainbuild 
+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
+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
+C
+C Add the components corresponding to local energy terms.
+C
+   10 continue
+      do i=1,nvar
+cd      write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
+        g(i)=g(i)+gloc(i,icg)
+      enddo
+C Uncomment following three lines for diagnostics.
+cd    call intout
+cd    call briefout(0,0.0d0)
+cd    write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n)
+      return
+      end
+C-------------------------------------------------------------------------
+      subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      include 'COMMON.IOUNITS'
+      external ufparm
+      integer uiparm(1)
+      double precision urparm(1)
+      dimension x(maxvar),g(maxvar)
+
+      icg=mod(nf,2)+1
+      if (nf-nfl+1) 20,30,40
+   20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm)
+c     write (iout,*) 'grad 20'
+      if (nf.eq.0) return
+      goto 40
+   30 continue
+#ifdef OSF
+c     Intercept NaNs in the coordinates
+c      write(iout,*) (var(i),i=1,nvar)
+      x_sum=0.D0
+      do i=1,n
+        x_sum=x_sum+x(i)
+      enddo
+      if (x_sum.ne.x_sum) then
+        write(iout,*)" *** grad_restr : Found NaN in coordinates"
+        call flush(iout)
+        print *," *** grad_restr : Found NaN in coordinates"
+        return
+      endif
+#endif
+      call var_to_geom_restr(n,x)
+      call chainbuild 
+C
+C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
+C
+   40 call cartder
+C
+C Convert the Cartesian gradient into internal-coordinate gradient.
+C
+
+      ig=0
+      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                                                                   
+        ig=ig+1
+        g(ig)=gphii
+       ELSE
+        ind=ind+nres-1-i
+       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
+       ENDIF
+      enddo
+
+      do i=2,nres-1
+       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
+         ENDIF
+        endif
+      enddo
+
+      
+      do i=2,nres-1
+        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
+         ENDIF
+        endif
+      enddo
+
+C
+C Add the components corresponding to local energy terms.
+C
+
+      ig=0
+      igall=0
+      do i=4,nres
+        igall=igall+1
+        if (mask_phi(i).eq.1) then
+          ig=ig+1
+          g(ig)=g(ig)+gloc(igall,icg)
+        endif
+      enddo
+
+      do i=3,nres
+        igall=igall+1
+        if (mask_theta(i).eq.1) then
+          ig=ig+1
+          g(ig)=g(ig)+gloc(igall,icg)
+        endif
+      enddo
+     
+      do ij=1,2
+      do i=2,nres-1
+        if (itype(i).ne.10) then
+          igall=igall+1
+          if (mask_side(i).eq.1) then
+            ig=ig+1
+            g(ig)=g(ig)+gloc(igall,icg)
+          endif
+        endif
+      enddo
+      enddo
+
+cd      do i=1,ig
+cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
+cd      enddo
+      return
+      end
+C-------------------------------------------------------------------------
+      subroutine cartgrad
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      include 'COMMON.MD'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.TIME1'
+      include 'COMMON.SCCOR'
+c
+c This subrouting calculates total Cartesian coordinate gradient. 
+c The subroutine chainbuild_cart and energy MUST be called beforehand.
+c
+c        do i=1,nres
+c        write (iout,*) "przed sum_grad", gloc_sc(1,i,icg),gloc(i,icg)
+c        enddo
+
+#ifdef TIMING
+      time00=MPI_Wtime()
+#endif
+      icg=1
+      call sum_gradient
+#ifdef TIMING
+#endif
+c        do i=1,nres
+c        write (iout,*) "checkgrad", gloc_sc(1,i,icg),gloc(i,icg)
+c        enddo     
+cd      write (iout,*) "After sum_gradient"
+cd      do i=1,nres-1
+cd        write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
+cd        write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
+cd      enddo
+c If performing constraint dynamics, add the gradients of the constraint energy
+      if(usampl.and.totT.gt.eq_time) then
+         do i=1,nct
+           do j=1,3
+             gradc(j,i,icg)=gradc(j,i,icg)+dudconst(j,i)+duscdiff(j,i)
+             gradx(j,i,icg)=gradx(j,i,icg)+dudxconst(j,i)+duscdiffx(j,i)
+           enddo
+         enddo
+         do i=1,nres-3
+           gloc(i,icg)=gloc(i,icg)+dugamma(i)
+         enddo
+         do i=1,nres-2
+           gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
+         enddo
+      endif 
+#ifdef TIMING
+      time01=MPI_Wtime()
+#endif
+      call intcartderiv
+#ifdef TIMING
+      time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
+#endif
+cd      call checkintcartgrad
+cd      write(iout,*) 'calling int_to_cart'
+cd      write (iout,*) "gcart, gxcart, gloc before int_to_cart"
+      do i=1,nct
+        do j=1,3
+          gcart(j,i)=gradc(j,i,icg)
+          gxcart(j,i)=gradx(j,i,icg)
+        enddo
+cd        write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3),
+cd     &    (gxcart(j,i),j=1,3),gloc(i,icg)
+      enddo
+#ifdef TIMING
+      time01=MPI_Wtime()
+#endif
+      call int_to_cart
+#ifdef TIMING
+      time_inttocart=time_inttocart+MPI_Wtime()-time01
+#endif
+cd      write (iout,*) "gcart and gxcart after int_to_cart"
+cd      do i=0,nres-1
+cd        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
+cd     &      (gxcart(j,i),j=1,3)
+cd      enddo
+#ifdef TIMING
+      time_cartgrad=time_cartgrad+MPI_Wtime()-time00
+#endif
+      return
+      end
+C-------------------------------------------------------------------------
+      subroutine zerograd
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.DERIV'
+      include 'COMMON.CHAIN'
+      include 'COMMON.VAR'
+      include 'COMMON.MD'
+      include 'COMMON.SCCOR'
+C
+C Initialize Cartesian-coordinate gradient
+C
+      do i=1,nres
+       do j=1,3
+         gvdwx(j,i)=0.0D0
+          gvdwxT(j,i)=0.0D0
+          gradx_scp(j,i)=0.0D0
+         gvdwc(j,i)=0.0D0
+          gvdwcT(j,i)=0.0D0
+          gvdwc_scp(j,i)=0.0D0
+          gvdwc_scpp(j,i)=0.0d0
+         gelc (j,i)=0.0D0
+         gelc_long(j,i)=0.0D0
+          gradb(j,i)=0.0d0
+          gradbx(j,i)=0.0d0
+          gvdwpp(j,i)=0.0d0
+          gel_loc(j,i)=0.0d0
+          gel_loc_long(j,i)=0.0d0
+         ghpbc(j,i)=0.0D0
+         ghpbx(j,i)=0.0D0
+          gcorr3_turn(j,i)=0.0d0
+          gcorr4_turn(j,i)=0.0d0
+          gradcorr(j,i)=0.0d0
+          gradcorr_long(j,i)=0.0d0
+          gradcorr5_long(j,i)=0.0d0
+          gradcorr6_long(j,i)=0.0d0
+          gcorr6_turn_long(j,i)=0.0d0
+          gradcorr5(j,i)=0.0d0
+          gradcorr6(j,i)=0.0d0
+          gcorr6_turn(j,i)=0.0d0
+          gsccorc(j,i)=0.0d0
+          gsccorx(j,i)=0.0d0
+          gradc(j,i,icg)=0.0d0
+          gradx(j,i,icg)=0.0d0
+          gscloc(j,i)=0.0d0
+          gsclocx(j,i)=0.0d0
+          do intertyp=1,3
+           gloc_sc(intertyp,i,icg)=0.0d0
+          enddo
+        enddo
+      enddo
+C
+C Initialize the gradient of local energy terms.
+C
+      do i=1,4*nres
+        gloc(i,icg)=0.0D0
+      enddo
+      do i=1,nres
+        gel_loc_loc(i)=0.0d0
+        gcorr_loc(i)=0.0d0
+        g_corr5_loc(i)=0.0d0
+        g_corr6_loc(i)=0.0d0
+        gel_loc_turn3(i)=0.0d0
+        gel_loc_turn4(i)=0.0d0
+        gel_loc_turn6(i)=0.0d0
+        gsccor_loc(i)=0.0d0
+      enddo
+c initialize gcart and gxcart
+      do i=0,nres
+        do j=1,3
+          gcart(j,i)=0.0d0
+          gxcart(j,i)=0.0d0
+        enddo
+      enddo
+      return
+      end
+c-------------------------------------------------------------------------
+      double precision function fdum()
+      fdum=0.0D0
+      return
+      end