source/unres/src-HCD-5D/check_ecartint_CASC_NC.F

   1       subroutine check_ecartint
   2 ! Check the gradient of the energy in Cartesian coordinates.
   3       implicit none
   4       include 'DIMENSIONS'
   5       include 'COMMON.CONTROL'
   6       include 'COMMON.CHAIN'
   7       include 'COMMON.INTERACT'
   8       include 'COMMON.DERIV'
   9       include 'COMMON.IOUNITS'
  10       include 'COMMON.VAR'
  11       include 'COMMON.CONTACTS'
  12       include 'COMMON.MD'
  13       include 'COMMON.LOCAL'
  14       include 'COMMON.SPLITELE'
  15       integer icall
  16       common /srutu/ icall
  17       double precision ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),
  18      & x(maxvar),g(maxvar)
  19       double precision dcnorm_safe(3),dxnorm_safe(3)
  20       double precision grad_s(6,0:maxres),grad_s1(6,0:maxres)
  21       double precision phi_temp(maxres),theta_temp(maxres),
  22      &  alph_temp(maxres),omeg_temp(maxres)
  23       double precision ddc1(3),ddcn(3),dcnorm_safe1(3),dcnorm_safe2(3)
  24       double precision energia(0:n_ene),energia1(0:n_ene)
  25       integer uiparm(1)
  26       double precision urparm(1)
  27       double precision fdum
  28       external fdum
  29       integer i,j,k,nf
  30       double precision etot,etot1,etot2,etot11,etot12,etot21,etot22
  31       double precision dist,alpha,beta
  32       icg=1
  33       nf=0
  34       nfl=0
  35       call intout
  36 !      call intcartderiv
  37 !      call checkintcartgrad
  38       call zerograd
  39       aincr=1.0D-5
  40       write(iout,*) 'Calling CHECK_ECARTINT.'
  41       nf=0
  42       icall=0
  43       write (iout,*) "Before geom_to_var"
  44       call geom_to_var(nvar,x)
  45       write (iout,*) "after geom_to_var"
  46       write (iout,*) "split_ene ",split_ene
  47       call flush(iout)
  48       if (.not.split_ene) then
  49         write(iout,*) 'Calling CHECK_ECARTINT if'
  50         call etotal(energia)
  51 !elwrite(iout,*) 'Calling CHECK_ECARTINT if'
  52         etot=energia(0)
  53         write (iout,*) "etot",etot
  54         call enerprint(energia(0))
  55         call flush(iout)
  56 !el        call enerprint(energia)
  57 !elwrite(iout,*) 'Calling CHECK_ECARTINT if'
  58 c        call flush(iout)
  59 c        write (iout,*) "enter cartgrad"
  60 c        call flush(iout)
  61         call cartgrad
  62 c Transform the gradient to the CA-SC basis
  63         call grad_transform
  64 !elwrite(iout,*) 'Calling CHECK_ECARTINT if'
  65 c        write (iout,*) "exit cartgrad"
  66 c        call flush(iout)
  67         icall =1
  68 c        do i=1,nres
  69 c          write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
  70 c        enddo
  71         do j=1,3
  72           grad_s(j,0)=gcart(j,0)
  73         enddo
  74 !elwrite(iout,*) 'Calling CHECK_ECARTINT if'
  75         do i=1,nres
  76           do j=1,3
  77             grad_s(j,i)=gcart(j,i)
  78             grad_s(j+3,i)=gxcart(j,i)
  79           enddo
  80         enddo
  81       else
  82 c        write(iout,*) 'Calling CHECK_ECARTIN else.'
  83 !- split gradient check
  84         call zerograd
  85         call etotal_long(energia)
  86         call enerprint(energia(0))
  87 !el        call enerprint(energia)
  88 c        call flush(iout)
  89 c        write (iout,*) "enter cartgrad"
  90 c        call flush(iout)
  91         call cartgrad
  92 c Transform the gradient to CA-SC coordinates
  93         call grad_transform
  94 c        write (iout,*) "exit cartgrad"
  95 c        call flush(iout)
  96         icall =1
  97         write (iout,*) "longrange grad"
  98         do i=1,nres
  99           write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 100      &    (gxcart(j,i),j=1,3)
 101         enddo
 102         do j=1,3
 103           grad_s(j,0)=gcart(j,0)
 104         enddo
 105         do i=1,nres
 106           do j=1,3
 107             grad_s(j,i)=gcart(j,i)
 108             grad_s(j+3,i)=gxcart(j,i)
 109           enddo
 110         enddo
 111         call zerograd
 112         call etotal_short(energia)
 113         call enerprint(energia(0))
 114         call flush(iout)
 115 c        write (iout,*) "enter cartgrad"
 116 c        call flush(iout)
 117         call cartgrad
 118 c        write (iout,*) "exit cartgrad"
 119 c        call flush(iout)
 120 c Transform the gradient to CA-SC basis
 121         call grad_transform
 122         icall =1
 123         write (iout,*) "shortrange grad"
 124         do i=1,nres
 125           write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 126      &    (gxcart(j,i),j=1,3)
 127         enddo
 128         do j=1,3
 129           grad_s1(j,0)=gcart(j,0)
 130         enddo
 131         do i=1,nres
 132           do j=1,3
 133             grad_s1(j,i)=gcart(j,i)
 134             grad_s1(j+3,i)=gxcart(j,i)
 135           enddo
 136         enddo
 137       endif
 138       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
 139 !      do i=1,nres
 140 c      do i=nnt,nct
 141       do i=1,nres
 142         do j=1,3
 143           if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1)
 144           if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres)
 145           ddc(j)=c(j,i)
 146           ddx(j)=c(j,i+nres)
 147           dcnorm_safe1(j)=dc_norm(j,i-1)
 148           dcnorm_safe2(j)=dc_norm(j,i)
 149           dxnorm_safe(j)=dc_norm(j,i+nres)
 150         enddo
 151         do j=1,3
 152           c(j,i)=ddc(j)+aincr
 153           if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=c(j,1)+aincr
 154           if (nct.lt.nres .and. i.eq.nct) c(j,nres)=c(j,nres)+aincr
 155           if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
 156           dc(j,i)=c(j,i+1)-c(j,i)
 157           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 158           call int_from_cart1(.false.)
 159           if (.not.split_ene) then
 160             call etotal(energia1)
 161             etot1=energia1(0)
 162 c            write (iout,*) "ij",i,j," etot1",etot1
 163           else
 164 !- split gradient
 165             call etotal_long(energia1)
 166             etot11=energia1(0)
 167             call etotal_short(energia1)
 168             etot12=energia1(0)
 169           endif
 170 !- end split gradient
 171 !          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
 172           c(j,i)=ddc(j)-aincr
 173           if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)-aincr
 174           if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)-aincr
 175           if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
 176           dc(j,i)=c(j,i+1)-c(j,i)
 177           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 178           call int_from_cart1(.false.)
 179           if (.not.split_ene) then
 180             call etotal(energia1)
 181             etot2=energia1(0)
 182 c            write (iout,*) "ij",i,j," etot2",etot2
 183             ggg(j)=(etot1-etot2)/(2*aincr)
 184           else
 185 !- split gradient
 186             call etotal_long(energia1)
 187             etot21=energia1(0)
 188             ggg(j)=(etot11-etot21)/(2*aincr)
 189             call etotal_short(energia1)
 190             etot22=energia1(0)
 191             ggg1(j)=(etot12-etot22)/(2*aincr)
 192 !- end split gradient
 193 !            write (iout,*) "etot21",etot21," etot22",etot22
 194           endif
 195 !          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
 196           c(j,i)=ddc(j)
 197           if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)
 198           if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)
 199           if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
 200           dc(j,i)=c(j,i+1)-c(j,i)
 201           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 202           dc_norm(j,i-1)=dcnorm_safe1(j)
 203           dc_norm(j,i)=dcnorm_safe2(j)
 204           dc_norm(j,i+nres)=dxnorm_safe(j)
 205         enddo
 206         do j=1,3
 207           c(j,i+nres)=ddx(j)+aincr
 208           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 209           call int_from_cart1(.false.)
 210           if (.not.split_ene) then
 211             call etotal(energia1)
 212             etot1=energia1(0)
 213           else
 214 !- split gradient
 215             call etotal_long(energia1)
 216             etot11=energia1(0)
 217             call etotal_short(energia1)
 218             etot12=energia1(0)
 219           endif
 220 !- end split gradient
 221           c(j,i+nres)=ddx(j)-aincr
 222           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 223           call int_from_cart1(.false.)
 224           if (.not.split_ene) then
 225             call etotal(energia1)
 226             etot2=energia1(0)
 227             ggg(j+3)=(etot1-etot2)/(2*aincr)
 228           else
 229 !- split gradient
 230             call etotal_long(energia1)
 231             etot21=energia1(0)
 232             ggg(j+3)=(etot11-etot21)/(2*aincr)
 233             call etotal_short(energia1)
 234             etot22=energia1(0)
 235             ggg1(j+3)=(etot12-etot22)/(2*aincr)
 236 !- end split gradient
 237           endif
 238 !          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
 239           c(j,i+nres)=ddx(j)
 240           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 241           dc_norm(j,i+nres)=dxnorm_safe(j)
 242           call int_from_cart1(.false.)
 243         enddo
 244         write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
 245      &   i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
 246         if (split_ene) then
 247           write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
 248      &   i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
 249      &   k=1,6)
 250          write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
 251      &   i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
 252      &   ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
 253         endif
 254       enddo
 255       return
 256       end