source/unres/src_MD-M/moments.f

   1       subroutine inertia_tensor
   2 c Calculating the intertia tensor for the entire protein in order to
   3 c remove the perpendicular components of velocity matrix which cause
   4 c the molecule to rotate.
   5        implicit real*8 (a-h,o-z)
   6        include 'DIMENSIONS'
   7        include 'COMMON.CONTROL'
   8        include 'COMMON.VAR'
   9        include 'COMMON.MD'
  10        include 'COMMON.CHAIN'
  11        include 'COMMON.DERIV'
  12        include 'COMMON.GEO'
  13        include 'COMMON.LOCAL'
  14        include 'COMMON.INTERACT'
  15        include 'COMMON.IOUNITS'
  16        include 'COMMON.NAMES'
  17
  18       double precision Im(3,3),Imcp(3,3),cm(3),pr(3),M_SC,
  19      & eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),
  20      & vpp(3,0:MAXRES),vs_p(3),pr1(3,3),
  21      & pr2(3,3),pp(3),incr(3),v(3),mag,mag2
  22       common /gucio/ cm
  23       integer iti,inres
  24         do i=1,3
  25            do j=1,3
  26              Im(i,j)=0.0d0
  27              pr1(i,j)=0.0d0
  28              pr2(i,j)=0.0d0
  29            enddo
  30            L(i)=0.0d0
  31            cm(i)=0.0d0
  32            vrot(i)=0.0d0
  33         enddo
  34 c   calculating the center of the mass of the protein
  35         do i=nnt,nct-1
  36           do j=1,3
  37             cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
  38           enddo
  39         enddo
  40         do j=1,3
  41          cm(j)=mp*cm(j)
  42         enddo
  43         M_SC=0.0d0
  44         do i=nnt,nct
  45            iti=iabs(itype(i))
  46            M_SC=M_SC+msc(iabs(iti))
  47            inres=i+nres
  48            do j=1,3
  49             cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
  50            enddo
  51         enddo
  52         do j=1,3
  53           cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
  54         enddo
  55
  56         do i=nnt,nct-1
  57           do j=1,3
  58             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
  59           enddo
  60           Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
  61           Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
  62           Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
  63           Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
  64           Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
  65           Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
  66         enddo
  67
  68         do i=nnt,nct
  69            iti=iabs(itype(i))
  70            inres=i+nres
  71            do j=1,3
  72              pr(j)=c(j,inres)-cm(j)
  73            enddo
  74           Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
  75           Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
  76           Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
  77           Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
  78           Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
  79           Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
  80         enddo
  81
  82         do i=nnt,nct-1
  83           Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))*
  84      &    vbld(i+1)*vbld(i+1)*0.25d0
  85           Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))*
  86      &    vbld(i+1)*vbld(i+1)*0.25d0
  87           Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))*
  88      &    vbld(i+1)*vbld(i+1)*0.25d0
  89           Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))*
  90      &    vbld(i+1)*vbld(i+1)*0.25d0
  91           Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))*
  92      &    vbld(i+1)*vbld(i+1)*0.25d0
  93           Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))*
  94      &    vbld(i+1)*vbld(i+1)*0.25d0
  95         enddo
  96
  97
  98         do i=nnt,nct
  99          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 100            iti=iabs(itype(i))
 101            inres=i+nres
 102           Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)*
 103      &    dc_norm(1,inres))*vbld(inres)*vbld(inres)
 104           Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)*
 105      &    dc_norm(2,inres))*vbld(inres)*vbld(inres)
 106           Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)*
 107      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 108           Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)*
 109      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 110           Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)*
 111      &    dc_norm(2,inres))*vbld(inres)*vbld(inres)
 112           Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)*
 113      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 114          endif
 115         enddo
 116
 117         call angmom(cm,L)
 118 c        write(iout,*) "The angular momentum before adjustment:"
 119 c        write(iout,*) (L(j),j=1,3)
 120
 121         Im(2,1)=Im(1,2)
 122         Im(3,1)=Im(1,3)
 123         Im(3,2)=Im(2,3)
 124
 125 c  Copying the Im matrix for the djacob subroutine
 126         do i=1,3
 127           do j=1,3
 128             Imcp(i,j)=Im(i,j)
 129             Id(i,j)=0.0d0
 130           enddo
 131         enddo
 132
 133 c   Finding the eigenvectors and eignvalues of the inertia tensor
 134 c       write (iout,*) "Calling djacob"
 135 c       call flush(iout)
 136        call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
 137 c       write (iout,*) "Eigenvalues & Eigenvectors"
 138 c       write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
 139 c       write (iout,*)
 140 c       do i=1,3
 141 c         write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
 142 c       enddo
 143 c       call flush(iout)
 144 c   Constructing the diagonalized matrix
 145        do i=1,3
 146          if (dabs(eigval(i)).gt.1.0d-15) then
 147            Id(i,i)=1.0d0/eigval(i)
 148          else
 149            Id(i,i)=0.0d0
 150          endif
 151        enddo
 152         do i=1,3
 153            do j=1,3
 154               Imcp(i,j)=eigvec(j,i)
 155            enddo
 156         enddo
 157         do i=1,3
 158            do j=1,3
 159               do k=1,3
 160                  pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
 161               enddo
 162            enddo
 163         enddo
 164         do i=1,3
 165            do j=1,3
 166               do k=1,3
 167                  pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
 168               enddo
 169            enddo
 170         enddo
 171 c  Calculating the total rotational velocity of the molecule
 172        do i=1,3
 173          do j=1,3
 174            vrot(i)=vrot(i)+pr2(i,j)*L(j)
 175          enddo
 176        enddo
 177 c   Resetting the velocities
 178        do i=nnt,nct-1
 179          call vecpr(vrot(1),dc(1,i),vp)
 180          do j=1,3
 181            d_t(j,i)=d_t(j,i)-vp(j)
 182           enddo
 183         enddo
 184         do i=nnt,nct
 185         if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 186            inres=i+nres
 187            call vecpr(vrot(1),dc(1,inres),vp)
 188            do j=1,3
 189              d_t(j,inres)=d_t(j,inres)-vp(j)
 190            enddo
 191         endif
 192        enddo
 193        call angmom(cm,L)
 194 c       write(iout,*) "The angular momentum after adjustment:"
 195 c       write(iout,*) (L(j),j=1,3)
 196        return
 197        end
 198 c----------------------------------------------------------------------------
 199        subroutine angmom(cm,L)
 200        implicit real*8 (a-h,o-z)
 201        include 'DIMENSIONS'
 202        include 'COMMON.CONTROL'
 203        include 'COMMON.VAR'
 204        include 'COMMON.MD'
 205        include 'COMMON.CHAIN'
 206        include 'COMMON.DERIV'
 207        include 'COMMON.GEO'
 208        include 'COMMON.LOCAL'
 209        include 'COMMON.INTERACT'
 210        include 'COMMON.IOUNITS'
 211        include 'COMMON.NAMES'
 212
 213       double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
 214      &  pp(3)
 215       integer iti,inres
 216 c  Calculate the angular momentum
 217        do j=1,3
 218           L(j)=0.0d0
 219        enddo
 220        do j=1,3
 221           incr(j)=d_t(j,0)
 222        enddo
 223        do i=nnt,nct-1
 224           do j=1,3
 225             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
 226           enddo
 227           do j=1,3
 228             v(j)=incr(j)+0.5d0*d_t(j,i)
 229           enddo
 230           do j=1,3
 231             incr(j)=incr(j)+d_t(j,i)
 232           enddo
 233           call vecpr(pr(1),v(1),vp)
 234           do j=1,3
 235             L(j)=L(j)+mp*vp(j)
 236           enddo
 237           do j=1,3
 238              pr(j)=0.5d0*dc(j,i)
 239              pp(j)=0.5d0*d_t(j,i)
 240           enddo
 241          call vecpr(pr(1),pp(1),vp)
 242          do j=1,3
 243              L(j)=L(j)+Ip*vp(j)
 244           enddo
 245         enddo
 246         do j=1,3
 247           incr(j)=d_t(j,0)
 248         enddo
 249         do i=nnt,nct
 250          iti=iabs(itype(i))
 251          inres=i+nres
 252          do j=1,3
 253            pr(j)=c(j,inres)-cm(j)
 254          enddo
 255          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 256            do j=1,3
 257              v(j)=incr(j)+d_t(j,inres)
 258            enddo
 259          else
 260            do j=1,3
 261              v(j)=incr(j)
 262            enddo
 263          endif
 264          call vecpr(pr(1),v(1),vp)
 265 c         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
 266 c     &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
 267          do j=1,3
 268             L(j)=L(j)+msc(iabs(iti))*vp(j)
 269          enddo
 270 c         write (iout,*) "L",(l(j),j=1,3)
 271          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 272            do j=1,3
 273             v(j)=incr(j)+d_t(j,inres)
 274            enddo
 275            call vecpr(dc(1,inres),d_t(1,inres),vp)
 276            do j=1,3
 277              L(j)=L(j)+Isc(iti)*vp(j)
 278           enddo
 279          endif
 280          do j=1,3
 281              incr(j)=incr(j)+d_t(j,i)
 282          enddo
 283        enddo
 284       return
 285       end
 286 c------------------------------------------------------------------------------
 287        subroutine vcm_vel(vcm)
 288        implicit real*8 (a-h,o-z)
 289        include 'DIMENSIONS'
 290        include 'COMMON.VAR'
 291        include 'COMMON.MD'
 292        include 'COMMON.CHAIN'
 293        include 'COMMON.DERIV'
 294        include 'COMMON.GEO'
 295        include 'COMMON.LOCAL'
 296        include 'COMMON.INTERACT'
 297        include 'COMMON.IOUNITS'
 298        double precision vcm(3),vv(3),summas,amas
 299        do j=1,3
 300          vcm(j)=0.0d0
 301          vv(j)=d_t(j,0)
 302        enddo
 303        summas=0.0d0
 304        do i=nnt,nct
 305          if (i.lt.nct) then
 306            summas=summas+mp
 307            do j=1,3
 308              vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
 309            enddo
 310          endif
 311          amas=msc(iabs(itype(i)))
 312          summas=summas+amas
 313          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 314            do j=1,3
 315              vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
 316            enddo
 317          else
 318            do j=1,3
 319              vcm(j)=vcm(j)+amas*vv(j)
 320            enddo
 321          endif
 322          do j=1,3
 323            vv(j)=vv(j)+d_t(j,i)
 324          enddo
 325        enddo
 326 c       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
 327        do j=1,3
 328          vcm(j)=vcm(j)/summas
 329        enddo
 330        return
 331        end