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 .or.
 187      &      itype(i).ne.ntyp1 .and. itype(i+1).eq.ntyp1) cycle
 188         call vecpr(vrot(1),dc(1,i),vp)
 189         do j=1,3
 190           d_t(j,i)=d_t(j,i)-vp(j)
 191         enddo
 192       enddo
 193 #else
 194       do i=nnt,nct-1
 195         call vecpr(vrot(1),dc(1,i),vp)
 196         do j=1,3
 197           d_t(j,i)=d_t(j,i)-vp(j)
 198         enddo
 199       enddo
 200 #endif
 201       do i=nnt,nct
 202         if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 203           inres=i+nres
 204           call vecpr(vrot(1),dc(1,inres),vp)
 205           do j=1,3
 206             d_t(j,inres)=d_t(j,inres)-vp(j)
 207           enddo
 208         endif
 209       enddo
 210       call angmom(cm,L)
 211 c       write(iout,*) "The angular momentum after adjustment:"
 212 c       write(iout,*) (L(j),j=1,3)
 213       return
 214       end
 215 c----------------------------------------------------------------------------
 216       subroutine angmom(cm,L)
 217       implicit none
 218       include 'DIMENSIONS'
 219       include 'COMMON.CONTROL'
 220       include 'COMMON.VAR'
 221       include 'COMMON.MD'
 222 #ifdef FIVEDIAG
 223       include 'COMMON.LAGRANGE.5diag'
 224 #else
 225       include 'COMMON.LAGRANGE'
 226 #endif
 227 #ifdef LANG0
 228 #ifdef FIVEDIAG
 229       include 'COMMON.LANGEVIN.lang0.5diag'
 230 #else
 231       include 'COMMON.LANGEVIN.lang0'
 232 #endif
 233 #else
 234       include 'COMMON.LANGEVIN'
 235 #endif
 236       include 'COMMON.CHAIN'
 237       include 'COMMON.DERIV'
 238       include 'COMMON.GEO'
 239       include 'COMMON.LOCAL'
 240       include 'COMMON.INTERACT'
 241       include 'COMMON.IOUNITS'
 242       include 'COMMON.NAMES'
 243       double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
 244      &  pp(3)
 245       integer iti,inres,i,j
 246 c  Calculate the angular momentum
 247       do j=1,3
 248          L(j)=0.0d0
 249       enddo
 250       do j=1,3
 251          incr(j)=d_t(j,0)
 252       enddo
 253       do i=nnt,nct-1
 254         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
 255         do j=1,3
 256           pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
 257         enddo
 258         do j=1,3
 259           v(j)=incr(j)+0.5d0*d_t(j,i)
 260         enddo
 261         do j=1,3
 262           incr(j)=incr(j)+d_t(j,i)
 263         enddo
 264         call vecpr(pr(1),v(1),vp)
 265         do j=1,3
 266           L(j)=L(j)+mp*vp(j)
 267         enddo
 268         do j=1,3
 269           pr(j)=0.5d0*dc(j,i)
 270           pp(j)=0.5d0*d_t(j,i)
 271         enddo
 272         call vecpr(pr(1),pp(1),vp)
 273         do j=1,3
 274           L(j)=L(j)+Ip*vp(j)
 275         enddo
 276       enddo
 277       do j=1,3
 278         incr(j)=d_t(j,0)
 279       enddo
 280       do i=nnt,nct
 281         iti=iabs(itype(i))
 282         inres=i+nres
 283         do j=1,3
 284           pr(j)=c(j,inres)-cm(j)
 285         enddo
 286         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 287           do j=1,3
 288             v(j)=incr(j)+d_t(j,inres)
 289           enddo
 290         else
 291           do j=1,3
 292             v(j)=incr(j)
 293           enddo
 294         endif
 295         call vecpr(pr(1),v(1),vp)
 296 c          write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
 297 c      &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
 298         do j=1,3
 299           L(j)=L(j)+msc(iabs(iti))*vp(j)
 300         enddo
 301 c          write (iout,*) "L",(l(j),j=1,3)
 302         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 303           do j=1,3
 304             v(j)=incr(j)+d_t(j,inres)
 305           enddo
 306           call vecpr(dc(1,inres),d_t(1,inres),vp)
 307           do j=1,3
 308             L(j)=L(j)+Isc(iti)*vp(j)
 309           enddo
 310         endif
 311         do j=1,3
 312             incr(j)=incr(j)+d_t(j,i)
 313         enddo
 314       enddo
 315       return
 316       end
 317 c------------------------------------------------------------------------------
 318        subroutine vcm_vel(vcm)
 319        implicit none
 320        include 'DIMENSIONS'
 321        include 'COMMON.VAR'
 322        include 'COMMON.MD'
 323 #ifdef FIVEDIAG
 324        include 'COMMON.LAGRANGE.5diag'
 325 #else
 326        include 'COMMON.LAGRANGE'
 327 #endif
 328        include 'COMMON.CHAIN'
 329        include 'COMMON.DERIV'
 330        include 'COMMON.GEO'
 331        include 'COMMON.LOCAL'
 332        include 'COMMON.INTERACT'
 333        include 'COMMON.IOUNITS'
 334        double precision vcm(3),vv(3),summas,amas
 335        integer i,j
 336        do j=1,3
 337          vcm(j)=0.0d0
 338          vv(j)=d_t(j,0)
 339        enddo
 340        summas=0.0d0
 341        do i=nnt,nct
 342          if (i.lt.nct) then
 343            summas=summas+mp
 344            do j=1,3
 345              vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
 346            enddo
 347          endif
 348          amas=msc(iabs(itype(i)))
 349          summas=summas+amas
 350          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 351            do j=1,3
 352              vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
 353            enddo
 354          else
 355            do j=1,3
 356              vcm(j)=vcm(j)+amas*vv(j)
 357            enddo
 358          endif
 359          do j=1,3
 360            vv(j)=vv(j)+d_t(j,i)
 361          enddo
 362        enddo
 363 c       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
 364        do j=1,3
 365          vcm(j)=vcm(j)/summas
 366        enddo
 367        return
 368        end
 369 #else
 370       subroutine inertia_tensor
 371 c Calculating the intertia tensor for the entire protein in order to
 372 c remove the perpendicular components of velocity matrix which cause
 373 c the molecule to rotate.
 374        implicit none
 375        include 'DIMENSIONS'
 376        include 'COMMON.CONTROL'
 377        include 'COMMON.VAR'
 378        include 'COMMON.MD'
 379 #ifdef FIVEDIAG
 380        include 'COMMON.LAGRANGE.5diag'
 381 #else
 382        include 'COMMON.LAGRANGE'
 383 #endif
 384        include 'COMMON.CHAIN'
 385        include 'COMMON.DERIV'
 386        include 'COMMON.GEO'
 387        include 'COMMON.LOCAL'
 388        include 'COMMON.INTERACT'
 389        include 'COMMON.IOUNITS'
 390        include 'COMMON.NAMES'
 391
 392       double precision Im(3,3),Imcp(3,3),cm(3),pr(3),M_SC,
 393      & eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),
 394      & vpp(3,0:MAXRES),vs_p(3),pr1(3,3),
 395      & pr2(3,3),pp(3),incr(3),v(3),mag,mag2
 396       common /gucio/ cm
 397       integer iti,inres,i,j,k
 398         do i=1,3
 399            do j=1,3
 400              Im(i,j)=0.0d0
 401              pr1(i,j)=0.0d0
 402              pr2(i,j)=0.0d0
 403            enddo
 404            L(i)=0.0d0
 405            cm(i)=0.0d0
 406            vrot(i)=0.0d0
 407         enddo
 408 c   calculating the center of the mass of the protein
 409         do i=nnt,nct-1
 410           do j=1,3
 411             cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
 412           enddo
 413         enddo
 414         do j=1,3
 415          cm(j)=mp*cm(j)
 416         enddo
 417         M_SC=0.0d0
 418         do i=nnt,nct
 419            iti=iabs(itype(i))
 420            M_SC=M_SC+msc(iabs(iti))
 421            inres=i+nres
 422            do j=1,3
 423             cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
 424            enddo
 425         enddo
 426         do j=1,3
 427           cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
 428         enddo
 429
 430         do i=nnt,nct-1
 431           do j=1,3
 432             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
 433           enddo
 434           Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
 435           Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
 436           Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
 437           Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
 438           Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
 439           Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
 440         enddo
 441
 442         do i=nnt,nct
 443            iti=iabs(itype(i))
 444            inres=i+nres
 445            do j=1,3
 446              pr(j)=c(j,inres)-cm(j)
 447            enddo
 448           Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
 449           Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
 450           Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
 451           Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
 452           Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
 453           Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
 454         enddo
 455
 456         do i=nnt,nct-1
 457           Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))*
 458      &    vbld(i+1)*vbld(i+1)*0.25d0
 459           Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))*
 460      &    vbld(i+1)*vbld(i+1)*0.25d0
 461           Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))*
 462      &    vbld(i+1)*vbld(i+1)*0.25d0
 463           Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))*
 464      &    vbld(i+1)*vbld(i+1)*0.25d0
 465           Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))*
 466      &    vbld(i+1)*vbld(i+1)*0.25d0
 467           Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))*
 468      &    vbld(i+1)*vbld(i+1)*0.25d0
 469         enddo
 470
 471
 472         do i=nnt,nct
 473          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 474            iti=iabs(itype(i))
 475            inres=i+nres
 476           Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)*
 477      &    dc_norm(1,inres))*vbld(inres)*vbld(inres)
 478           Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)*
 479      &    dc_norm(2,inres))*vbld(inres)*vbld(inres)
 480           Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)*
 481      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 482           Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)*
 483      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 484           Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)*
 485      &    dc_norm(2,inres))*vbld(inres)*vbld(inres)
 486           Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)*
 487      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
 488          endif
 489         enddo
 490
 491         call angmom(cm,L)
 492 c        write(iout,*) "The angular momentum before adjustment:"
 493 c        write(iout,*) (L(j),j=1,3)
 494
 495         Im(2,1)=Im(1,2)
 496         Im(3,1)=Im(1,3)
 497         Im(3,2)=Im(2,3)
 498
 499 c  Copying the Im matrix for the djacob subroutine
 500         do i=1,3
 501           do j=1,3
 502             Imcp(i,j)=Im(i,j)
 503             Id(i,j)=0.0d0
 504           enddo
 505         enddo
 506
 507 c   Finding the eigenvectors and eignvalues of the inertia tensor
 508        call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
 509 c       write (iout,*) "Eigenvalues & Eigenvectors"
 510 c       write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
 511 c       write (iout,*)
 512 c       do i=1,3
 513 c         write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
 514 c       enddo
 515 c   Constructing the diagonalized matrix
 516        do i=1,3
 517          if (dabs(eigval(i)).gt.1.0d-15) then
 518            Id(i,i)=1.0d0/eigval(i)
 519          else
 520            Id(i,i)=0.0d0
 521          endif
 522        enddo
 523         do i=1,3
 524            do j=1,3
 525               Imcp(i,j)=eigvec(j,i)
 526            enddo
 527         enddo
 528         do i=1,3
 529            do j=1,3
 530               do k=1,3
 531                  pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
 532               enddo
 533            enddo
 534         enddo
 535         do i=1,3
 536            do j=1,3
 537               do k=1,3
 538                  pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
 539               enddo
 540            enddo
 541         enddo
 542 c  Calculating the total rotational velocity of the molecule
 543        do i=1,3
 544          do j=1,3
 545            vrot(i)=vrot(i)+pr2(i,j)*L(j)
 546          enddo
 547        enddo
 548 c   Resetting the velocities
 549        do i=nnt,nct-1
 550          call vecpr(vrot(1),dc(1,i),vp)
 551          do j=1,3
 552            d_t(j,i)=d_t(j,i)-vp(j)
 553           enddo
 554         enddo
 555         do i=nnt,nct
 556         if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 557            inres=i+nres
 558            call vecpr(vrot(1),dc(1,inres),vp)
 559            do j=1,3
 560              d_t(j,inres)=d_t(j,inres)-vp(j)
 561            enddo
 562         endif
 563        enddo
 564        call angmom(cm,L)
 565 c       write(iout,*) "The angular momentum after adjustment:"
 566 c       write(iout,*) (L(j),j=1,3)
 567        return
 568        end
 569 c----------------------------------------------------------------------------
 570        subroutine angmom(cm,L)
 571        implicit none
 572        include 'DIMENSIONS'
 573        include 'COMMON.CONTROL'
 574        include 'COMMON.VAR'
 575        include 'COMMON.MD'
 576 #ifdef FIVEDIAG
 577        include 'COMMON.LAGRANGE.5diag'
 578 #else
 579        include 'COMMON.LAGRANGE'
 580 #endif
 581 #ifdef LANG0
 582 #ifdef FIVEDIAG
 583       include 'COMMON.LANGEVIN.lang0.5diag'
 584 #else
 585       include 'COMMON.LANGEVIN.lang0'
 586 #endif
 587 #else
 588        include 'COMMON.LANGEVIN'
 589 #endif
 590        include 'COMMON.CHAIN'
 591        include 'COMMON.DERIV'
 592        include 'COMMON.GEO'
 593        include 'COMMON.LOCAL'
 594        include 'COMMON.INTERACT'
 595        include 'COMMON.IOUNITS'
 596        include 'COMMON.NAMES'
 597
 598       double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
 599      &  pp(3)
 600       integer iti,inres,i,j
 601 c  Calculate the angular momentum
 602        do j=1,3
 603           L(j)=0.0d0
 604        enddo
 605        do j=1,3
 606           incr(j)=d_t(j,0)
 607        enddo
 608        do i=nnt,nct-1
 609           do j=1,3
 610             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
 611           enddo
 612           do j=1,3
 613             v(j)=incr(j)+0.5d0*d_t(j,i)
 614           enddo
 615           do j=1,3
 616             incr(j)=incr(j)+d_t(j,i)
 617           enddo
 618           call vecpr(pr(1),v(1),vp)
 619           do j=1,3
 620             L(j)=L(j)+mp*vp(j)
 621           enddo
 622           do j=1,3
 623              pr(j)=0.5d0*dc(j,i)
 624              pp(j)=0.5d0*d_t(j,i)
 625           enddo
 626          call vecpr(pr(1),pp(1),vp)
 627          do j=1,3
 628              L(j)=L(j)+Ip*vp(j)
 629           enddo
 630         enddo
 631         do j=1,3
 632           incr(j)=d_t(j,0)
 633         enddo
 634         do i=nnt,nct
 635          iti=iabs(itype(i))
 636          inres=i+nres
 637          do j=1,3
 638            pr(j)=c(j,inres)-cm(j)
 639          enddo
 640          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 641            do j=1,3
 642              v(j)=incr(j)+d_t(j,inres)
 643            enddo
 644          else
 645            do j=1,3
 646              v(j)=incr(j)
 647            enddo
 648          endif
 649          call vecpr(pr(1),v(1),vp)
 650 c         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
 651 c     &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
 652          do j=1,3
 653             L(j)=L(j)+msc(iabs(iti))*vp(j)
 654          enddo
 655 c         write (iout,*) "L",(l(j),j=1,3)
 656          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 657            do j=1,3
 658             v(j)=incr(j)+d_t(j,inres)
 659            enddo
 660            call vecpr(dc(1,inres),d_t(1,inres),vp)
 661            do j=1,3
 662              L(j)=L(j)+Isc(iti)*vp(j)
 663           enddo
 664          endif
 665          do j=1,3
 666              incr(j)=incr(j)+d_t(j,i)
 667          enddo
 668        enddo
 669       return
 670       end
 671 c------------------------------------------------------------------------------
 672        subroutine vcm_vel(vcm)
 673        implicit none
 674        include 'DIMENSIONS'
 675        include 'COMMON.VAR'
 676        include 'COMMON.MD'
 677 #ifdef FIVEDIAG
 678        include 'COMMON.LAGRANGE.5diag'
 679 #else
 680        include 'COMMON.LAGRANGE'
 681 #endif
 682        include 'COMMON.CHAIN'
 683        include 'COMMON.DERIV'
 684        include 'COMMON.GEO'
 685        include 'COMMON.LOCAL'
 686        include 'COMMON.INTERACT'
 687        include 'COMMON.IOUNITS'
 688        double precision vcm(3),vv(3),summas,amas
 689        integer i,j
 690        do j=1,3
 691          vcm(j)=0.0d0
 692          vv(j)=d_t(j,0)
 693        enddo
 694        summas=0.0d0
 695        do i=nnt,nct
 696          if (i.lt.nct) then
 697            summas=summas+mp
 698            do j=1,3
 699              vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
 700            enddo
 701          endif
 702          amas=msc(iabs(itype(i)))
 703          summas=summas+amas
 704          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 705            do j=1,3
 706              vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
 707            enddo
 708          else
 709            do j=1,3
 710              vcm(j)=vcm(j)+amas*vv(j)
 711            enddo
 712          endif
 713          do j=1,3
 714            vv(j)=vv(j)+d_t(j,i)
 715          enddo
 716        enddo
 717 c       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
 718        do j=1,3
 719          vcm(j)=vcm(j)/summas
 720        enddo
 721        return
 722        end
 723 #endif