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        call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
 135 c       write (iout,*) "Eigenvalues & Eigenvectors"
 136 c       write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
 137 c       write (iout,*)
 138 c       do i=1,3
 139 c         write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
 140 c       enddo
 141 c   Constructing the diagonalized matrix
 142        do i=1,3
 143          if (dabs(eigval(i)).gt.1.0d-15) then
 144            Id(i,i)=1.0d0/eigval(i)
 145          else
 146            Id(i,i)=0.0d0
 147          endif
 148        enddo
 149         do i=1,3
 150            do j=1,3
 151               Imcp(i,j)=eigvec(j,i)
 152            enddo
 153         enddo
 154         do i=1,3
 155            do j=1,3
 156               do k=1,3
 157                  pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
 158               enddo
 159            enddo
 160         enddo
 161         do i=1,3
 162            do j=1,3
 163               do k=1,3
 164                  pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
 165               enddo
 166            enddo
 167         enddo
 168 c  Calculating the total rotational velocity of the molecule
 169        do i=1,3
 170          do j=1,3
 171            vrot(i)=vrot(i)+pr2(i,j)*L(j)
 172          enddo
 173        enddo
 174 c   Resetting the velocities
 175        do i=nnt,nct-1
 176          call vecpr(vrot(1),dc(1,i),vp)
 177          do j=1,3
 178            d_t(j,i)=d_t(j,i)-vp(j)
 179           enddo
 180         enddo
 181         do i=nnt,nct
 182         if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 183            inres=i+nres
 184            call vecpr(vrot(1),dc(1,inres),vp)
 185            do j=1,3
 186              d_t(j,inres)=d_t(j,inres)-vp(j)
 187            enddo
 188         endif
 189        enddo
 190        call angmom(cm,L)
 191 c       write(iout,*) "The angular momentum after adjustment:"
 192 c       write(iout,*) (L(j),j=1,3)
 193        return
 194        end
 195 c----------------------------------------------------------------------------
 196        subroutine angmom(cm,L)
 197        implicit real*8 (a-h,o-z)
 198        include 'DIMENSIONS'
 199        include 'COMMON.CONTROL'
 200        include 'COMMON.VAR'
 201        include 'COMMON.MD'
 202        include 'COMMON.CHAIN'
 203        include 'COMMON.DERIV'
 204        include 'COMMON.GEO'
 205        include 'COMMON.LOCAL'
 206        include 'COMMON.INTERACT'
 207        include 'COMMON.IOUNITS'
 208        include 'COMMON.NAMES'
 209
 210       double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
 211      &  pp(3)
 212       integer iti,inres
 213 c  Calculate the angular momentum
 214        do j=1,3
 215           L(j)=0.0d0
 216        enddo
 217        do j=1,3
 218           incr(j)=d_t(j,0)
 219        enddo
 220        do i=nnt,nct-1
 221           do j=1,3
 222             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
 223           enddo
 224           do j=1,3
 225             v(j)=incr(j)+0.5d0*d_t(j,i)
 226           enddo
 227           do j=1,3
 228             incr(j)=incr(j)+d_t(j,i)
 229           enddo
 230           call vecpr(pr(1),v(1),vp)
 231           do j=1,3
 232             L(j)=L(j)+mp*vp(j)
 233           enddo
 234           do j=1,3
 235              pr(j)=0.5d0*dc(j,i)
 236              pp(j)=0.5d0*d_t(j,i)
 237           enddo
 238          call vecpr(pr(1),pp(1),vp)
 239          do j=1,3
 240              L(j)=L(j)+Ip*vp(j)
 241           enddo
 242         enddo
 243         do j=1,3
 244           incr(j)=d_t(j,0)
 245         enddo
 246         do i=nnt,nct
 247          iti=iabs(itype(i))
 248          inres=i+nres
 249          do j=1,3
 250            pr(j)=c(j,inres)-cm(j)
 251          enddo
 252          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 253            do j=1,3
 254              v(j)=incr(j)+d_t(j,inres)
 255            enddo
 256          else
 257            do j=1,3
 258              v(j)=incr(j)
 259            enddo
 260          endif
 261          call vecpr(pr(1),v(1),vp)
 262 c         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
 263 c     &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
 264          do j=1,3
 265             L(j)=L(j)+msc(iabs(iti))*vp(j)
 266          enddo
 267 c         write (iout,*) "L",(l(j),j=1,3)
 268          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 269            do j=1,3
 270             v(j)=incr(j)+d_t(j,inres)
 271            enddo
 272            call vecpr(dc(1,inres),d_t(1,inres),vp)
 273            do j=1,3
 274              L(j)=L(j)+Isc(iti)*vp(j)
 275           enddo
 276          endif
 277          do j=1,3
 278              incr(j)=incr(j)+d_t(j,i)
 279          enddo
 280        enddo
 281       return
 282       end
 283 c------------------------------------------------------------------------------
 284        subroutine vcm_vel(vcm)
 285        implicit real*8 (a-h,o-z)
 286        include 'DIMENSIONS'
 287        include 'COMMON.VAR'
 288        include 'COMMON.MD'
 289        include 'COMMON.CHAIN'
 290        include 'COMMON.DERIV'
 291        include 'COMMON.GEO'
 292        include 'COMMON.LOCAL'
 293        include 'COMMON.INTERACT'
 294        include 'COMMON.IOUNITS'
 295        double precision vcm(3),vv(3),summas,amas
 296        do j=1,3
 297          vcm(j)=0.0d0
 298          vv(j)=d_t(j,0)
 299        enddo
 300        summas=0.0d0
 301        do i=nnt,nct
 302          if (i.lt.nct) then
 303            summas=summas+mp
 304            do j=1,3
 305              vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
 306            enddo
 307          endif
 308          amas=msc(iabs(itype(i)))
 309          summas=summas+amas
 310          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 311            do j=1,3
 312              vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
 313            enddo
 314          else
 315            do j=1,3
 316              vcm(j)=vcm(j)+amas*vv(j)
 317            enddo
 318          endif
 319          do j=1,3
 320            vv(j)=vv(j)+d_t(j,i)
 321          enddo
 322        enddo
 323 c       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
 324        do j=1,3
 325          vcm(j)=vcm(j)/summas
 326        enddo
 327        return
 328        end