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

   1 #ifdef FIVEDIAG
   2       subroutine inertia_tensor
   3 c Calculating the intertia tensor for the entire protein in order to
   4 c remove the perpendicular components of velocity matrix which cause
   5 c the molecule to rotate.
   6        implicit none
   7        include 'DIMENSIONS'
   8        include 'COMMON.CONTROL'
   9        include 'COMMON.VAR'
  10        include 'COMMON.MD'
  11 #ifdef FIVEDIAG
  12        include 'COMMON.LAGRANGE.5diag'
  13 #else
  14        include 'COMMON.LAGRANGE'
  15 #endif
  16        include 'COMMON.CHAIN'
  17        include 'COMMON.DERIV'
  18        include 'COMMON.GEO'
  19        include 'COMMON.LOCAL'
  20        include 'COMMON.INTERACT'
  21        include 'COMMON.IOUNITS'
  22        include 'COMMON.NAMES'
  23
  24       double precision Im(3,3),Imcp(3,3),cm(3),pr(3),M_SC,
  25      & eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),
  26      & vpp(3,0:MAXRES),vs_p(3),pr1(3,3),
  27      & pr2(3,3),pp(3),incr(3),v(3),mag,mag2
  28       common /gucio/ cm
  29       integer iti,inres,i,j,k
  30       do i=1,3
  31         do j=1,3
  32           Im(i,j)=0.0d0
  33           pr1(i,j)=0.0d0
  34           pr2(i,j)=0.0d0
  35         enddo
  36         L(i)=0.0d0
  37         cm(i)=0.0d0
  38         vrot(i)=0.0d0
  39       enddo
  40 c   caulating the center of the mass of the protein
  41       do i=nnt,nct-1
  42         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
  43         do j=1,3
  44           cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
  45         enddo
  46       enddo
  47       do j=1,3
  48        cm(j)=mp*cm(j)
  49       enddo
  50       M_SC=0.0d0
  51       do i=nnt,nct
  52          iti=iabs(itype(i))
  53          if (iti.eq.ntyp1) cycle
  54          M_SC=M_SC+msc(iabs(iti))
  55          inres=i+nres
  56          do j=1,3
  57           cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
  58          enddo
  59       enddo
  60       do j=1,3
  61         cm(j)=cm(j)/(M_SC+dimenp*mp)
  62       enddo
  63
  64       do i=nnt,nct-1
  65         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
  66         do j=1,3
  67           pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
  68         enddo
  69         Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
  70         Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
  71         Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
  72         Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
  73         Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
  74         Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
  75       enddo
  76
  77       do i=nnt,nct
  78         iti=iabs(itype(i))
  79         if (iti.eq.ntyp1) cycle
  80         inres=i+nres
  81         do j=1,3
  82           pr(j)=c(j,inres)-cm(j)
  83         enddo
  84         Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
  85         Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
  86         Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
  87         Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
  88         Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
  89         Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
  90       enddo
  91
  92       do i=nnt,nct-1
  93         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
  94         Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))*
  95      &  vbld(i+1)*vbld(i+1)*0.25d0
  96         Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))*
  97      &  vbld(i+1)*vbld(i+1)*0.25d0
  98         Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))*
  99      &  vbld(i+1)*vbld(i+1)*0.25d0
 100         Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))*
 101      &  vbld(i+1)*vbld(i+1)*0.25d0
 102         Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))*
 103      &  vbld(i+1)*vbld(i+1)*0.25d0
 104         Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))*
 105      &  vbld(i+1)*vbld(i+1)*0.25d0
 106       enddo
 107       do i=nnt,nct
 108         iti=iabs(itype(i))
 109         if (iti.ne.10 .and. iti.ne.ntyp1) then
 110           inres=i+nres
 111           Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)*
 112      &    dc_norm(1,inres))*vbld(inres)*vbld(inres)
 113           Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)*
 114      &    dc_norm(2,inres))*vbld(inres)*vbld(inres)
 115           Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)*
 116      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 117           Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)*
 118      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 119           Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)*
 120      &    dc_norm(2,inres))*vbld(inres)*vbld(inres)
 121           Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)*
 122      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 123         endif
 124       enddo
 125
 126       call angmom(cm,L)
 127 c      write(iout,*) "The angular momentum before adjustment:"
 128 c      write(iout,*) (L(j),j=1,3)
 129
 130       Im(2,1)=Im(1,2)
 131       Im(3,1)=Im(1,3)
 132       Im(3,2)=Im(2,3)
 133
 134 c  Copng the Im matrix for the djacob subroutine
 135       do i=1,3
 136         do j=1,3
 137           Imcp(i,j)=Im(i,j)
 138           Id(i,j)=0.0d0
 139         enddo
 140       enddo
 141 c   Finding the eigenvectors and eignvalues of the inertia tensor
 142       call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
 143 c       write (iout,*) "Eigenvalues & Eigenvectors"
 144 c       write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
 145 c       write (iout,*)
 146 c       do i=1,3
 147 c         write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
 148 c       enddo
 149 c   Constructing the diagonalized matrix
 150       do i=1,3
 151         if (dabs(eigval(i)).gt.1.0d-15) then
 152           Id(i,i)=1.0d0/eigval(i)
 153         else
 154           Id(i,i)=0.0d0
 155         endif
 156       enddo
 157       do i=1,3
 158         do j=1,3
 159           Imcp(i,j)=eigvec(j,i)
 160         enddo
 161       enddo
 162       do i=1,3
 163         do j=1,3
 164           do k=1,3
 165              pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
 166           enddo
 167         enddo
 168       enddo
 169       do i=1,3
 170         do j=1,3
 171           do k=1,3
 172             pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
 173           enddo
 174         enddo
 175       enddo
 176 c  Calculating the total rotational velocity of the molecule
 177       do i=1,3
 178         do j=1,3
 179           vrot(i)=vrot(i)+pr2(i,j)*L(j)
 180         enddo
 181       enddo
 182 c   Resetting the velocities
 183 #ifdef FIVEDIAG
 184       do i=nnt,nct-1
 185         write (iout,*) itype(i+1),itype(i)
 186         if (itype(i+1).ne.ntyp1 .and. itype(i).eq.ntyp1) cycle
 187         call vecpr(vrot(1),dc(1,i),vp)
 188         do j=1,3
 189           d_t(j,i)=d_t(j,i)-vp(j)
 190         enddo
 191       enddo
 192 #else
 193       do i=nnt,nct-1
 194         call vecpr(vrot(1),dc(1,i),vp)
 195         do j=1,3
 196           d_t(j,i)=d_t(j,i)-vp(j)
 197         enddo
 198       enddo
 199 #endif
 200       do i=nnt,nct
 201         if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 202           inres=i+nres
 203           call vecpr(vrot(1),dc(1,inres),vp)
 204           do j=1,3
 205             d_t(j,inres)=d_t(j,inres)-vp(j)
 206           enddo
 207         endif
 208       enddo
 209       call angmom(cm,L)
 210 c       write(iout,*) "The angular momentum after adjustment:"
 211 c       write(iout,*) (L(j),j=1,3)
 212       return
 213       end
 214 c----------------------------------------------------------------------------
 215       subroutine angmom(cm,L)
 216       implicit none
 217       include 'DIMENSIONS'
 218       include 'COMMON.CONTROL'
 219       include 'COMMON.VAR'
 220       include 'COMMON.MD'
 221 #ifdef FIVEDIAG
 222       include 'COMMON.LAGRANGE.5diag'
 223 #else
 224       include 'COMMON.LAGRANGE'
 225 #endif
 226 #ifdef LANG0
 227 #ifdef FIVEDIAG
 228       include 'COMMON.LANGEVIN.lang0.5diag'
 229 #else
 230       include 'COMMON.LANGEVIN.lang0'
 231 #endif
 232 #else
 233       include 'COMMON.LANGEVIN'
 234 #endif
 235       include 'COMMON.CHAIN'
 236       include 'COMMON.DERIV'
 237       include 'COMMON.GEO'
 238       include 'COMMON.LOCAL'
 239       include 'COMMON.INTERACT'
 240       include 'COMMON.IOUNITS'
 241       include 'COMMON.NAMES'
 242       double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
 243      &  pp(3)
 244       integer iti,inres,i,j
 245 c  Calculate the angular momentum
 246       do j=1,3
 247          L(j)=0.0d0
 248       enddo
 249       do j=1,3
 250          incr(j)=d_t(j,0)
 251       enddo
 252       do i=nnt,nct-1
 253         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
 254         do j=1,3
 255           pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
 256         enddo
 257         do j=1,3
 258           v(j)=incr(j)+0.5d0*d_t(j,i)
 259         enddo
 260         do j=1,3
 261           incr(j)=incr(j)+d_t(j,i)
 262         enddo
 263         call vecpr(pr(1),v(1),vp)
 264         do j=1,3
 265           L(j)=L(j)+mp*vp(j)
 266         enddo
 267         do j=1,3
 268           pr(j)=0.5d0*dc(j,i)
 269           pp(j)=0.5d0*d_t(j,i)
 270         enddo
 271         call vecpr(pr(1),pp(1),vp)
 272         do j=1,3
 273           L(j)=L(j)+Ip*vp(j)
 274         enddo
 275       enddo
 276       do j=1,3
 277         incr(j)=d_t(j,0)
 278       enddo
 279       do i=nnt,nct
 280         iti=iabs(itype(i))
 281         inres=i+nres
 282         do j=1,3
 283           pr(j)=c(j,inres)-cm(j)
 284         enddo
 285         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 286           do j=1,3
 287             v(j)=incr(j)+d_t(j,inres)
 288           enddo
 289         else
 290           do j=1,3
 291             v(j)=incr(j)
 292           enddo
 293         endif
 294         call vecpr(pr(1),v(1),vp)
 295 c          write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
 296 c      &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
 297         do j=1,3
 298           L(j)=L(j)+msc(iabs(iti))*vp(j)
 299         enddo
 300 c          write (iout,*) "L",(l(j),j=1,3)
 301         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 302           do j=1,3
 303             v(j)=incr(j)+d_t(j,inres)
 304           enddo
 305           call vecpr(dc(1,inres),d_t(1,inres),vp)
 306           do j=1,3
 307             L(j)=L(j)+Isc(iti)*vp(j)
 308           enddo
 309         endif
 310         do j=1,3
 311             incr(j)=incr(j)+d_t(j,i)
 312         enddo
 313       enddo
 314       return
 315       end
 316 c------------------------------------------------------------------------------
 317        subroutine vcm_vel(vcm)
 318        implicit none
 319        include 'DIMENSIONS'
 320        include 'COMMON.VAR'
 321        include 'COMMON.MD'
 322 #ifdef FIVEDIAG
 323        include 'COMMON.LAGRANGE.5diag'
 324 #else
 325        include 'COMMON.LAGRANGE'
 326 #endif
 327        include 'COMMON.CHAIN'
 328        include 'COMMON.DERIV'
 329        include 'COMMON.GEO'
 330        include 'COMMON.LOCAL'
 331        include 'COMMON.INTERACT'
 332        include 'COMMON.IOUNITS'
 333        double precision vcm(3),vv(3),summas,amas
 334        integer i,j
 335        do j=1,3
 336          vcm(j)=0.0d0
 337          vv(j)=d_t(j,0)
 338        enddo
 339        summas=0.0d0
 340        do i=nnt,nct
 341          if (i.lt.nct) then
 342            summas=summas+mp
 343            do j=1,3
 344              vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
 345            enddo
 346          endif
 347          amas=msc(iabs(itype(i)))
 348          summas=summas+amas
 349          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 350            do j=1,3
 351              vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
 352            enddo
 353          else
 354            do j=1,3
 355              vcm(j)=vcm(j)+amas*vv(j)
 356            enddo
 357          endif
 358          do j=1,3
 359            vv(j)=vv(j)+d_t(j,i)
 360          enddo
 361        enddo
 362 c       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
 363        do j=1,3
 364          vcm(j)=vcm(j)/summas
 365        enddo
 366        return
 367        end
 368 #else
 369       subroutine inertia_tensor
 370 c Calculating the intertia tensor for the entire protein in order to
 371 c remove the perpendicular components of velocity matrix which cause
 372 c the molecule to rotate.
 373        implicit none
 374        include 'DIMENSIONS'
 375        include 'COMMON.CONTROL'
 376        include 'COMMON.VAR'
 377        include 'COMMON.MD'
 378 #ifdef FIVEDIAG
 379        include 'COMMON.LAGRANGE.5diag'
 380 #else
 381        include 'COMMON.LAGRANGE'
 382 #endif
 383        include 'COMMON.CHAIN'
 384        include 'COMMON.DERIV'
 385        include 'COMMON.GEO'
 386        include 'COMMON.LOCAL'
 387        include 'COMMON.INTERACT'
 388        include 'COMMON.IOUNITS'
 389        include 'COMMON.NAMES'
 390
 391       double precision Im(3,3),Imcp(3,3),cm(3),pr(3),M_SC,
 392      & eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),
 393      & vpp(3,0:MAXRES),vs_p(3),pr1(3,3),
 394      & pr2(3,3),pp(3),incr(3),v(3),mag,mag2
 395       common /gucio/ cm
 396       integer iti,inres,i,j,k
 397         do i=1,3
 398            do j=1,3
 399              Im(i,j)=0.0d0
 400              pr1(i,j)=0.0d0
 401              pr2(i,j)=0.0d0
 402            enddo
 403            L(i)=0.0d0
 404            cm(i)=0.0d0
 405            vrot(i)=0.0d0
 406         enddo
 407 c   calculating the center of the mass of the protein
 408         do i=nnt,nct-1
 409           do j=1,3
 410             cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
 411           enddo
 412         enddo
 413         do j=1,3
 414          cm(j)=mp*cm(j)
 415         enddo
 416         M_SC=0.0d0
 417         do i=nnt,nct
 418            iti=iabs(itype(i))
 419            M_SC=M_SC+msc(iabs(iti))
 420            inres=i+nres
 421            do j=1,3
 422             cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
 423            enddo
 424         enddo
 425         do j=1,3
 426           cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
 427         enddo
 428
 429         do i=nnt,nct-1
 430           do j=1,3
 431             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
 432           enddo
 433           Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
 434           Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
 435           Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
 436           Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
 437           Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
 438           Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
 439         enddo
 440
 441         do i=nnt,nct
 442            iti=iabs(itype(i))
 443            inres=i+nres
 444            do j=1,3
 445              pr(j)=c(j,inres)-cm(j)
 446            enddo
 447           Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
 448           Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
 449           Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
 450           Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
 451           Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
 452           Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
 453         enddo
 454
 455         do i=nnt,nct-1
 456           Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))*
 457      &    vbld(i+1)*vbld(i+1)*0.25d0
 458           Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))*
 459      &    vbld(i+1)*vbld(i+1)*0.25d0
 460           Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))*
 461      &    vbld(i+1)*vbld(i+1)*0.25d0
 462           Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))*
 463      &    vbld(i+1)*vbld(i+1)*0.25d0
 464           Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))*
 465      &    vbld(i+1)*vbld(i+1)*0.25d0
 466           Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))*
 467      &    vbld(i+1)*vbld(i+1)*0.25d0
 468         enddo
 469
 470
 471         do i=nnt,nct
 472          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 473            iti=iabs(itype(i))
 474            inres=i+nres
 475           Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)*
 476      &    dc_norm(1,inres))*vbld(inres)*vbld(inres)
 477           Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)*
 478      &    dc_norm(2,inres))*vbld(inres)*vbld(inres)
 479           Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)*
 480      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 481           Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)*
 482      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 483           Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)*
 484      &    dc_norm(2,inres))*vbld(inres)*vbld(inres)
 485           Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)*
 486      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 487          endif
 488         enddo
 489
 490         call angmom(cm,L)
 491 c        write(iout,*) "The angular momentum before adjustment:"
 492 c        write(iout,*) (L(j),j=1,3)
 493
 494         Im(2,1)=Im(1,2)
 495         Im(3,1)=Im(1,3)
 496         Im(3,2)=Im(2,3)
 497
 498 c  Copying the Im matrix for the djacob subroutine
 499         do i=1,3
 500           do j=1,3
 501             Imcp(i,j)=Im(i,j)
 502             Id(i,j)=0.0d0
 503           enddo
 504         enddo
 505
 506 c   Finding the eigenvectors and eignvalues of the inertia tensor
 507        call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
 508 c       write (iout,*) "Eigenvalues & Eigenvectors"
 509 c       write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
 510 c       write (iout,*)
 511 c       do i=1,3
 512 c         write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
 513 c       enddo
 514 c   Constructing the diagonalized matrix
 515        do i=1,3
 516          if (dabs(eigval(i)).gt.1.0d-15) then
 517            Id(i,i)=1.0d0/eigval(i)
 518          else
 519            Id(i,i)=0.0d0
 520          endif
 521        enddo
 522         do i=1,3
 523            do j=1,3
 524               Imcp(i,j)=eigvec(j,i)
 525            enddo
 526         enddo
 527         do i=1,3
 528            do j=1,3
 529               do k=1,3
 530                  pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
 531               enddo
 532            enddo
 533         enddo
 534         do i=1,3
 535            do j=1,3
 536               do k=1,3
 537                  pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
 538               enddo
 539            enddo
 540         enddo
 541 c  Calculating the total rotational velocity of the molecule
 542        do i=1,3
 543          do j=1,3
 544            vrot(i)=vrot(i)+pr2(i,j)*L(j)
 545          enddo
 546        enddo
 547 c   Resetting the velocities
 548        do i=nnt,nct-1
 549          call vecpr(vrot(1),dc(1,i),vp)
 550          do j=1,3
 551            d_t(j,i)=d_t(j,i)-vp(j)
 552           enddo
 553         enddo
 554         do i=nnt,nct
 555         if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 556            inres=i+nres
 557            call vecpr(vrot(1),dc(1,inres),vp)
 558            do j=1,3
 559              d_t(j,inres)=d_t(j,inres)-vp(j)
 560            enddo
 561         endif
 562        enddo
 563        call angmom(cm,L)
 564 c       write(iout,*) "The angular momentum after adjustment:"
 565 c       write(iout,*) (L(j),j=1,3)
 566        return
 567        end
 568 c----------------------------------------------------------------------------
 569        subroutine angmom(cm,L)
 570        implicit none
 571        include 'DIMENSIONS'
 572        include 'COMMON.CONTROL'
 573        include 'COMMON.VAR'
 574        include 'COMMON.MD'
 575 #ifdef FIVEDIAG
 576        include 'COMMON.LAGRANGE.5diag'
 577 #else
 578        include 'COMMON.LAGRANGE'
 579 #endif
 580 #ifdef LANG0
 581 #ifdef FIVEDIAG
 582       include 'COMMON.LANGEVIN.lang0.5diag'
 583 #else
 584       include 'COMMON.LANGEVIN.lang0'
 585 #endif
 586 #else
 587        include 'COMMON.LANGEVIN'
 588 #endif
 589        include 'COMMON.CHAIN'
 590        include 'COMMON.DERIV'
 591        include 'COMMON.GEO'
 592        include 'COMMON.LOCAL'
 593        include 'COMMON.INTERACT'
 594        include 'COMMON.IOUNITS'
 595        include 'COMMON.NAMES'
 596
 597       double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
 598      &  pp(3)
 599       integer iti,inres,i,j
 600 c  Calculate the angular momentum
 601        do j=1,3
 602           L(j)=0.0d0
 603        enddo
 604        do j=1,3
 605           incr(j)=d_t(j,0)
 606        enddo
 607        do i=nnt,nct-1
 608           do j=1,3
 609             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
 610           enddo
 611           do j=1,3
 612             v(j)=incr(j)+0.5d0*d_t(j,i)
 613           enddo
 614           do j=1,3
 615             incr(j)=incr(j)+d_t(j,i)
 616           enddo
 617           call vecpr(pr(1),v(1),vp)
 618           do j=1,3
 619             L(j)=L(j)+mp*vp(j)
 620           enddo
 621           do j=1,3
 622              pr(j)=0.5d0*dc(j,i)
 623              pp(j)=0.5d0*d_t(j,i)
 624           enddo
 625          call vecpr(pr(1),pp(1),vp)
 626          do j=1,3
 627              L(j)=L(j)+Ip*vp(j)
 628           enddo
 629         enddo
 630         do j=1,3
 631           incr(j)=d_t(j,0)
 632         enddo
 633         do i=nnt,nct
 634          iti=iabs(itype(i))
 635          inres=i+nres
 636          do j=1,3
 637            pr(j)=c(j,inres)-cm(j)
 638          enddo
 639          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 640            do j=1,3
 641              v(j)=incr(j)+d_t(j,inres)
 642            enddo
 643          else
 644            do j=1,3
 645              v(j)=incr(j)
 646            enddo
 647          endif
 648          call vecpr(pr(1),v(1),vp)
 649 c         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
 650 c     &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
 651          do j=1,3
 652             L(j)=L(j)+msc(iabs(iti))*vp(j)
 653          enddo
 654 c         write (iout,*) "L",(l(j),j=1,3)
 655          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 656            do j=1,3
 657             v(j)=incr(j)+d_t(j,inres)
 658            enddo
 659            call vecpr(dc(1,inres),d_t(1,inres),vp)
 660            do j=1,3
 661              L(j)=L(j)+Isc(iti)*vp(j)
 662           enddo
 663          endif
 664          do j=1,3
 665              incr(j)=incr(j)+d_t(j,i)
 666          enddo
 667        enddo
 668       return
 669       end
 670 c------------------------------------------------------------------------------
 671        subroutine vcm_vel(vcm)
 672        implicit none
 673        include 'DIMENSIONS'
 674        include 'COMMON.VAR'
 675        include 'COMMON.MD'
 676 #ifdef FIVEDIAG
 677        include 'COMMON.LAGRANGE.5diag'
 678 #else
 679        include 'COMMON.LAGRANGE'
 680 #endif
 681        include 'COMMON.CHAIN'
 682        include 'COMMON.DERIV'
 683        include 'COMMON.GEO'
 684        include 'COMMON.LOCAL'
 685        include 'COMMON.INTERACT'
 686        include 'COMMON.IOUNITS'
 687        double precision vcm(3),vv(3),summas,amas
 688        integer i,j
 689        do j=1,3
 690          vcm(j)=0.0d0
 691          vv(j)=d_t(j,0)
 692        enddo
 693        summas=0.0d0
 694        do i=nnt,nct
 695          if (i.lt.nct) then
 696            summas=summas+mp
 697            do j=1,3
 698              vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
 699            enddo
 700          endif
 701          amas=msc(iabs(itype(i)))
 702          summas=summas+amas
 703          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 704            do j=1,3
 705              vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
 706            enddo
 707          else
 708            do j=1,3
 709              vcm(j)=vcm(j)+amas*vv(j)
 710            enddo
 711          endif
 712          do j=1,3
 713            vv(j)=vv(j)+d_t(j,i)
 714          enddo
 715        enddo
 716 c       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
 717        do j=1,3
 718          vcm(j)=vcm(j)/summas
 719        enddo
 720        return
 721        end
 722 #endif