source/unres/REMD.F90

   1       module REMD
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use MD_data
   5       use REMD_data
   6       use muca_md
   7
   8       implicit none
   9 !-----------------------------------------------------------------------------
  10 !
  11 !
  12 !-----------------------------------------------------------------------------
  13       contains
  14 !-----------------------------------------------------------------------------
  15 ! lagrangian_lesyng.F
  16 !-----------------------------------------------------------------------------
  17       subroutine lagrangian
  18 !-------------------------------------------------------------------------
  19 !  This subroutine contains the total lagrangain from which the accelerations
  20 !  are obtained.  For numerical gradient checking, the derivetive of the
  21 !  lagrangian in the velocities and coordinates are calculated seperately
  22 !-------------------------------------------------------------------------
  23 !       implicit real*8 (a-h,o-z)
  24 !       include 'DIMENSIONS'
  25       use comm_cipiszcze
  26       use energy_data
  27 #ifdef FIVEDIAG
  28       use energy, only: grad_transform
  29 #endif
  30       use geometry_data, only: nres
  31       use control_data    !el, only: mucadyn,lmuca
  32 #ifdef MPI
  33        include 'mpif.h'
  34       real(kind=8) :: time00
  35 #endif
  36 !       include 'COMMON.VAR'
  37 !       include 'COMMON.CHAIN'
  38 !       include 'COMMON.DERIV'
  39 !       include 'COMMON.GEO'
  40 !       include 'COMMON.LOCAL'
  41 !       include 'COMMON.INTERACT'
  42 !       include 'COMMON.MD'
  43 !       include 'COMMON.IOUNITS'
  44 !       include 'COMMON.CONTROL'
  45 !       include 'COMMON.MUCA'
  46 !      include 'COMMON.TIME1'
  47        integer ::i,j,ind,itime,mnum,innt,inct,inct_prev,ichain,n,mark
  48        real(kind=8) :: zapas(6*nres) !,muca_factor      !maxres6=6*maxres
  49 ! the line below might be wrong
  50 #ifdef FIVEDIAG
  51        real(kind=8) :: rs(2*nres),xsolv(2*nres)
  52 #ifdef CHECK5DSOL
  53        real(kind=8) :: rscheck(2*nres),rsold(2*nres)
  54 #endif
  55 #endif
  56        logical :: lprn = .false.
  57 !el       common /cipiszcze/ itime
  58        itime = itt_comm
  59
  60 #ifdef TIMING
  61        time00=MPI_Wtime()
  62 #endif
  63 #ifdef FIVEDIAG
  64       call grad_transform
  65       d_a=0.0d0
  66       if (lprn) then
  67         write (iout,*) "Potential forces backbone"
  68         do i=1,nres
  69           write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3)
  70         enddo
  71         write (iout,*) "Potential forces sidechain"
  72         do i=nnt,nct
  73 !          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
  74           write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
  75         enddo
  76       endif
  77       do ichain=1,nchain
  78         n=dimen_chain(ichain)
  79         innt=iposd_chain(ichain)
  80         do j=1,3
  81           ind=1
  82           do i=chain_border(1,ichain),chain_border(2,ichain)
  83             mnum=molnum(i)
  84             if (itype(i,1).eq.10.or.mnum.ge.3)then
  85               rs(ind)=-gcart(j,i)-gxcart(j,i)
  86               ind=ind+1
  87             else
  88               rs(ind)=-gcart(j,i)
  89               rs(ind+1)=-gxcart(j,i)
  90               ind=ind+2
  91             end if
  92           enddo
  93
  94 #ifdef CHECK5DSOL
  95          rsold=rs
  96 #endif
  97          if (lprn) then
  98            write(iout,*) "RHS of the 5-diag equations system",&
  99            ichain," j",j
 100            do i=1,n
 101              write(iout,*) i,rs(i)
 102            enddo
 103          endif
 104
 105          call FDISYS (n,DM(innt),DU1(innt),DU2(innt),rs,xsolv)
 106          if (lprn) then
 107            write (iout,*) "Solution of the 5-diagonal equations system"
 108            do i=1,n
 109              write (iout,'(i5,f10.5)') i,xsolv(i)
 110            enddo
 111          endif
 112 #ifdef CHECK5DSOL
 113 ! Check the solution
 114           call fivediagmult(n,DMorig(innt),DU1orig(innt),DU2orig(innt),
 115      &      xsolv,rscheck)
 116           do i=1,n
 117             write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),
 118      &       "ratio",rscheck(i)/rsold(i)
 119           enddo
 120 ! end check
 121 #endif
 122 #undef CHECK5DSOL
 123           ind=1
 124           do i=chain_border(1,ichain),chain_border(2,ichain)
 125             mnum=molnum(i)
 126             if (itype(i,1).eq.10 .or.mnum.ge.3) then
 127               d_a(j,i)=xsolv(ind)
 128               ind=ind+1
 129             else
 130               d_a(j,i)=xsolv(ind)
 131               d_a(j,i+nres)=xsolv(ind+1)
 132               ind=ind+2
 133             end if
 134           enddo
 135         enddo ! j
 136       enddo ! ichain
 137       if (lprn) then
 138         write (iout,*) "Acceleration in CA and SC oordinates"
 139         do i=1,nres
 140           write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
 141         enddo
 142         do i=nnt,nct
 143           write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
 144         enddo
 145       endif
 146 !C Conevert d_a to virtual-bon-vector basis
 147 #define WLOS
 148 #ifdef WLOS
 149 !c      write (iout,*) "WLOS"
 150       if (nnt.eq.1) then
 151         d_a(:,0)=d_a(:,1)
 152       endif
 153       do i=1,nres
 154         mnum=molnum(i)
 155         if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum) .or.mnum.ge.3) then
 156           do j=1,3
 157             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 158           enddo
 159         else
 160           do j=1,3
 161             d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
 162             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 163           enddo
 164         end if
 165       enddo
 166       d_a(:,nres)=0.0d0
 167       d_a(:,nct)=0.0d0
 168       d_a(:,2*nres)=0.0d0
 169 !c      d_a(:,0)=d_a(:,1)
 170 !c      d_a(:,1)=0.0d0
 171 !c      write (iout,*) "Shifting accelerations"
 172       if (nnt.gt.1) then
 173         d_a(:,0)=d_a(:,1)
 174         d_a(:,1)=0.0d0
 175       endif
 176 #define CHUJ
 177 #ifdef CHUJ
 178       do ichain=2,nchain
 179 !c        write (iout,*) "ichain",chain_border1(1,ichain)-1,
 180 !c     &     chain_border1(1,ichain)
 181         d_a(:,chain_border1(1,ichain)-1)=d_a(:,chain_border1(1,ichain))
 182         d_a(:,chain_border1(1,ichain))=0.0d0
 183       enddo
 184 !c      write (iout,*) "Adding accelerations"
 185       do ichain=2,nchain
 186 !c        write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
 187 !c     &   chain_border(2,ichain-1)
 188         d_a(:,chain_border1(1,ichain)-1)=&
 189        d_a(:,chain_border1(1,ichain)-1)+d_a(:,chain_border(2,ichain-1))
 190         d_a(:,chain_border(2,ichain-1))=0.0d0
 191       enddo
 192 #endif
 193 #else
 194       inct_prev=0
 195       do j=1,3
 196         aaux(j)=0.0d0
 197       enddo
 198       do ichain=1,nchain
 199         innt=chain_border(1,ichain)
 200         inct=chain_border(2,ichain)
 201         do j=1,3
 202           d_a(j,inct_prev)=d_a(j,innt)-aaux(j)
 203         enddo
 204         inct_prev=inct+1
 205         do i=innt,inct
 206           mnum=molnum(i)
 207           if ((itype(i).ne.10).and.(mnum.le.3)) then
 208             do j=1,3
 209               d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
 210             enddo
 211           endif
 212         enddo
 213         do j=1,3
 214           aaux(j)=d_a(j,inct)
 215         enddo
 216         do i=innt,inct
 217           do j=1,3
 218             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 219           enddo
 220         enddo
 221       enddo
 222 #endif
 223       if (lprn) then
 224         write(iout,*) 'acceleration 3D FIVEDIAG in dC and dX'
 225         do i=0,nres
 226           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 227         enddo
 228         do i=nnt,nct
 229           write (iout,'(i3,3f10.5,3x,3f10.5)')&
 230           i,(d_a(j,i+nres),j=1,3)
 231         enddo
 232       endif
 233
 234 #else
 235 ! Old procedure
 236        do j=1,3
 237          zapas(j)=-gcart(j,0)
 238        enddo
 239       ind=3
 240       if (lprn) then
 241         write (iout,*) "Potential forces backbone"
 242       endif
 243       do i=nnt,nct-1
 244         if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
 245           i,(-gcart(j,i),j=1,3)
 246         do j=1,3
 247           ind=ind+1
 248           zapas(ind)=-gcart(j,i)
 249         enddo
 250       enddo
 251       if (lprn) write (iout,*) "Potential forces sidechain"
 252       do i=nnt,nct
 253           mnum=molnum(i)
 254         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
 255           .and.mnum.lt.4) then
 256           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
 257              i,(-gxcart(j,i),j=1,3)
 258           do j=1,3
 259             ind=ind+1
 260             zapas(ind)=-gxcart(j,i)
 261           enddo
 262         endif
 263       enddo
 264
 265       call ginv_mult(zapas,d_a_work)
 266
 267       do j=1,3
 268         d_a(j,0)=d_a_work(j)
 269       enddo
 270       ind=3
 271       do i=nnt,nct-1
 272         do j=1,3
 273           ind=ind+1
 274           d_a(j,i)=d_a_work(ind)
 275         enddo
 276       enddo
 277       do i=nnt,nct
 278        mnum=molnum(i)
 279 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
 280          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
 281           .and.mnum.lt.4) then
 282           do j=1,3
 283             ind=ind+1
 284             d_a(j,i+nres)=d_a_work(ind)
 285           enddo
 286         endif
 287       enddo
 288
 289 #endif
 290       if(lmuca) then
 291        imtime=imtime+1
 292        if(mucadyn.gt.0) call muca_update(potE)
 293        factor=muca_factor(potE)*t_bath*Rb
 294
 295 !d       print *,'lmuca ',factor,potE
 296        do j=1,3
 297           d_a(j,0)=d_a(j,0)*factor
 298        enddo
 299        do i=nnt,nct-1
 300          do j=1,3
 301           d_a(j,i)=d_a(j,i)*factor
 302          enddo
 303        enddo
 304        do i=nnt,nct
 305          do j=1,3
 306           d_a(j,i+nres)=d_a(j,i+nres)*factor
 307          enddo
 308        enddo
 309
 310       endif
 311
 312       if (lprn) then
 313         write(iout,*) 'acceleration 3D'
 314         write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
 315         do i=nnt,nct-1
 316          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 317         enddo
 318         do i=nnt,nct
 319          write (iout,'(i3,3f10.5,3x,3f10.5)') &
 320            i+nres,(d_a(j,i+nres),j=1,3)
 321         enddo
 322       endif
 323 #ifdef TIMING
 324       time_lagrangian=time_lagrangian+MPI_Wtime()-time00
 325 #endif
 326       return
 327       end subroutine lagrangian
 328 !-----------------------------------------------------------------------------
 329       subroutine setup_MD_matrices
 330
 331       use geometry_data, only: nres,nside
 332       use control_data
 333       use MPI_data
 334       use energy_data
 335       use geometry, only:int_bounds
 336       use md_calc
 337 !      implicit real*8 (a-h,o-z)
 338 !      include 'DIMENSIONS'
 339 #ifdef MPI
 340       include 'mpif.h'
 341       integer :: ierror
 342       real(kind=8) :: time00
 343 !      real(kind=8) ,allocatable, dimension(:)  :: DDM,DDU1,DDU2
 344 #endif
 345 !      include 'COMMON.SETUP'
 346 !      include 'COMMON.VAR'
 347 !      include 'COMMON.CHAIN'
 348 !      include 'COMMON.DERIV'
 349 !      include 'COMMON.GEO'
 350 !      include 'COMMON.LOCAL'
 351 !      include 'COMMON.INTERACT'
 352 !      include 'COMMON.MD'
 353 !#ifndef LANG0
 354 !      include 'COMMON.LANGEVIN'
 355 !#else
 356 !      include 'COMMON.LANGEVIN.lang0'
 357 !#endif
 358 !      include 'COMMON.IOUNITS'
 359 !      include 'COMMON.TIME1'
 360       logical :: lprn = .false.
 361       logical :: osob
 362 #ifndef FIVEDIAG
 363       real(kind=8),allocatable,dimension(:,:) :: Bmat,matmult
 364 #endif
 365       real(kind=8) :: dtdi
 366       real(kind=8),dimension(:),allocatable :: massvec,sqreig   !(maxres2) maxres2=2*maxres
 367       real(kind=8) :: relfeh,eps1,eps2
 368 !el      real(kind=8),dimension(:),allocatable :: Ghalf
 369 !el      real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf   !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
 370 !el      real(kind=8),dimension(2*nres,2*nres) :: Gcopy !(maxres2,maxres2)
 371 !el      real(kind=8),dimension(:,:),allocatable :: Gcopy
 372       real(kind=8),dimension(:),allocatable :: work     !(8*maxres6)
 373       integer,dimension(:),allocatable :: iwork !(maxres6) maxres6=6*maxres
 374 !      common /jakistam/ iwork,work,massvec,sqreig
 375 !el      common /przechowalnia/ Gcopy,Ghalf
 376       real(kind=8) :: coeff,mscab
 377       integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
 378       real(kind=8) :: ip4
 379       integer :: iz,mnum,ichain,n,dimenp,innt,inct
 380       if(.not.allocated(massvec)) then
 381       allocate(massvec(2*nres),sqreig(2*nres))
 382       allocate(work(8*6*nres))
 383       allocate(iwork(6*nres))
 384       endif
 385       print *,"just entered"
 386       relfeh=1.0d-14
 387       nres2=2*nres
 388       print *,"FIVEDIAG"
 389       write (iout,*) "before FIVEDIAG"
 390 #ifdef FIVEDIAG
 391       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
 392 #endif
 393 #ifndef FIVEDIAG
 394       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
 395       write (iout,*) "ALLOCATE"
 396       print *,"ALLOCATE"
 397       if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
 398 !      if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
 399       if(.not.allocated(Bmat))  allocate(Bmat(nres2,nres2))
 400       if(.not.allocated(matmult)) allocate(matmult(nres2,nres2))
 401 #endif
 402 !
 403 ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
 404 ! inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
 405 !
 406 ! Determine the number of degrees of freedom (dimen) and the number of
 407 ! sites (dimen1)
 408 !      dimen=(nct-nnt+1)+nside
 409 !      dimen1=(nct-nnt)+(nct-nnt+1)
 410 !      dimen3=dimen*3
 411 !      write (iout,*) "nnt",nnt," nct",nct," nside",nside
 412 #ifdef FIVEDIAG
 413 #ifdef CRYST_BOND
 414        ip4=ip/4
 415        DM(1)=mp/4+ip4
 416        if (iabs(itype(nnt)).eq.10) then
 417          DM(1)=DM(1)+msc(10)
 418          ind=2
 419          nind=1
 420        else
 421          DM(1)=DM(1)+isc(iabs(itype(nnt)))
 422          DM(2)=msc(iabs(itype(nnt)))+isc(iabs(itype(nnt)))
 423          ind=3
 424          nind=2
 425        endif
 426        do i=nnt+1,nct-1
 427         mnum=molnum(i)
 428 !         if (iabs(itype(i,1)).eq.ntyp1) cycle
 429          DM(ind)=2*ip4+mp/2
 430         if (iabs(itype(i,1)).eq.10 .or. &
 431           iabs(itype(i,mnum)).eq.ntyp1_molec(mnum) .or. mnum.ge.4) then
 432            if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10)
 433            ind=ind+1
 434          else
 435            DM(ind)=DM(ind)+isc(iabs(itype(i,1)))
 436            DM(ind+1)=msc(iabs(itype(i,1)))+isc(iabs(itype(i,1)))
 437            ind=ind+2
 438          endif
 439        enddo
 440        if (nct.gt.nnt) then
 441        DM(ind)=ip4+mp/4
 442        if (iabs(itype(nct)).eq.10) then
 443          DM(ind)=DM(ind)+msc(10)
 444          nind=ind
 445        else
 446          DM(ind)=DM(ind)+isc(iabs(itype(nct)))
 447          DM(ind+1)=msc(iabs(itype(nct)))+isc(iabs(itype(nct)))
 448          nind=ind+1
 449         endif
 450         endif
 451
 452
 453         ind=1
 454         do i=nnt,nct
 455         mnum=molnum(i)
 456         if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 &
 457           .and.mnum.lt.4) then
 458           DU1(ind)=-isc(iabs(itype(i,1)))
 459           DU1(ind+1)=0.0d0
 460           ind=ind+2
 461         else
 462           DU1(ind)=mp/4-ip4
 463           ind=ind+1
 464         endif
 465         enddo
 466
 467        ind=1
 468        do i=nnt,nct-1
 469        mnum=molnum(i)
 470 !        if (iabs(itype(i,1)).eq.ntyp1) cycle
 471         write (iout,*) "i",i," itype",itype(i,1),ntyp1
 472         if (iabs(itype(i,1)).ne.10 .and. &
 473           iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.lt.4) then
 474          DU2(ind)=mp/4-ip4
 475          DU2(ind+1)=0.0d0
 476          ind=ind+2
 477         else
 478          DU2(ind)=0.0d0
 479          DU2(ind+1)=0.0d0
 480          ind=ind+1
 481         endif
 482        enddo
 483 #else
 484       dimen=0
 485       dimen1=0
 486       dimenp=0
 487       do ichain=1,nchain
 488         dimen=dimen+chain_length(ichain)
 489         dimen1=dimen1+2*chain_length(ichain)-1
 490         dimenp=dimenp+chain_length(ichain)-1
 491       enddo
 492       write (iout,*) "Number of Calphas",dimen
 493       write (iout,*) "Number of sidechains",nside
 494       write (iout,*) "Number of peptide groups",dimenp
 495       dimen=dimen+nside ! number of centers
 496       dimen3=3*dimen  ! degrees of freedom
 497       write (iout,*) "Number of centers",dimen
 498       write (iout,*) "Degrees of freedom:",dimen3
 499 !      ip4=ip/4
 500       ind=1
 501       do ichain=1,nchain
 502         iposd_chain(ichain)=ind
 503         innt=chain_border(1,ichain)
 504         mnum=molnum(innt)
 505         inct=chain_border(2,ichain)
 506         if (mnum.eq.5) mp(mnum)=0.0
 507 !        if (mnum.eq.5) mp(mnum)=msc(itype(innt,mnum),mnum)
 508         DM(ind)=mp(mnum)/4+ip(mnum)/4
 509         if (iabs(itype(innt,1)).eq.10.or.molnum(innt).gt.2) then
 510            DM(ind)=DM(ind)+msc(itype(innt,mnum),mnum)
 511           ind=ind+1
 512           nind=1
 513         else
 514           DM(ind)=DM(ind)+isc(iabs(itype(innt,mnum)),mnum)
 515           DM(ind+1)=msc(iabs(itype(innt,mnum)),mnum)+isc(iabs(itype(innt,mnum)),mnum)
 516           ind=ind+2
 517           nind=2
 518         endif
 519         write (iout,*) "ind",ind," nind",nind
 520         do i=innt+1,inct-1
 521          mnum=molnum(i)
 522 !        if (iabs(itype(i)).eq.ntyp1) cycle
 523 !          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 524 !          if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum)
 525         if (mnum.eq.5) mp(mnum)=0.0
 526 !        if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 527           DM(ind)=2*ip(mnum)/4+mp(mnum)/2
 528           if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2) then
 529 !            if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2)&
 530 !             DM(ind)=DM(ind)+msc(itype(i,molnum(i)),mnum)
 531              DM(ind)=DM(ind)+msc(itype(i,mnum),mnum)
 532             ind=ind+1
 533             nind=nind+1
 534           else
 535             DM(ind)=DM(ind)+isc(iabs(itype(i,mnum)),mnum)
 536             DM(ind+1)=msc(iabs(itype(i,mnum)),mnum)+isc(iabs(itype(i,mnum)),mnum)
 537             ind=ind+2
 538             nind=nind+2
 539           endif
 540           write (iout,*) "i",i," ind",ind," nind",nind
 541         enddo
 542         if (inct.gt.innt) then
 543 !         DM(ind)=ip4+mp(molnum(inct))/4
 544           mnum=molnum(inct)
 545         if (mnum.eq.5) mp(mnum)=0.0
 546 !        if (mnum.eq.5) mp(mnum)=msc(itype(inct,molnum(inct)),molnum(inct))
 547
 548           DM(ind)=mp(mnum)/4+ip(mnum)/4
 549           if (iabs(itype(inct,mnum)).eq.10.or.molnum(inct).gt.2) then
 550               DM(ind)=DM(ind)+msc(itype(inct,molnum(inct)),molnum(inct))
 551             ind=ind+1
 552             nind=nind+1
 553           else
 554             mnum=molnum(inct)
 555             DM(ind)=DM(ind)+isc(iabs(itype(inct,mnum)),mnum)
 556             DM(ind+1)=msc(iabs(itype(inct,mnum)),mnum)+isc(iabs(itype(inct,mnum)),mnum)
 557             ind=ind+2
 558             nind=nind+2
 559           endif
 560         endif
 561         write (iout,*) "ind",ind," nind",nind
 562         dimen_chain(ichain)=nind
 563         enddo
 564
 565       do ichain=1,nchain
 566         ind=iposd_chain(ichain)
 567         innt=chain_border(1,ichain)
 568         inct=chain_border(2,ichain)
 569         do i=innt,inct
 570          mnum=molnum(i)
 571           if (iabs(itype(i,1)).ne.10 .and.iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.le.2) then
 572             DU1(ind)=-isc(iabs(itype(i,mnum)),mnum)
 573             DU1(ind+1)=0.0d0
 574             ind=ind+2
 575           else
 576 !           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 577 !            if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum)
 578              if (mnum.eq.5) mp(mnum)=0.0
 579             DU1(ind)=mp(mnum)/4-ip(mnum)/4
 580             ind=ind+1
 581           endif
 582         enddo
 583       enddo
 584
 585       do ichain=1,nchain
 586         ind=iposd_chain(ichain)
 587         innt=chain_border(1,ichain)
 588         inct=chain_border(2,ichain)
 589         do i=innt,inct-1
 590          mnum=molnum(i)
 591 !       if (iabs(itype(i)).eq.ntyp1) cycle
 592 !c          write (iout,*) "i",i," itype",itype(i),ntyp1
 593           if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.le.2) then
 594 !            if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 595         if (mnum.eq.5) mp(mnum)=0.0
 596             DU2(ind)=mp(mnum)/4-ip(mnum)/4
 597             DU2(ind+1)=0.0d0
 598             ind=ind+2
 599           else
 600             DU2(ind)=0.0d0
 601             DU2(ind+1)=0.0d0
 602             ind=ind+1
 603           endif
 604         enddo
 605       enddo
 606       DMorig=DM
 607       DU1orig=DU1
 608       DU2orig=DU2
 609       if (gmatout) then
 610       write (iout,*)"The upper part of the five-diagonal inertia matrix"
 611       endif
 612       do ichain=1,nchain
 613         if (gmatout) write (iout,'(a,i5)') 'Chain',ichain
 614         n=dimen_chain(ichain)
 615         innt=iposd_chain(ichain)
 616         inct=iposd_chain(ichain)+dimen_chain(ichain)-1
 617         if (gmatout) then
 618         do i=innt,inct
 619           if (i.lt.inct-1) then
 620             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
 621           else if (i.eq.inct-1) then
 622             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
 623           else
 624             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
 625           endif
 626         enddo
 627         endif
 628         call FDISYP (n, DM(innt:inct), DU1(innt:inct-1),&
 629         DU2(innt:inct-1), MARK)
 630
 631         if (mark.eq.-1) then
 632           write(iout,*)&
 633         "ERROR: the inertia matrix is not positive definite for chain",&
 634          ichain
 635 #ifdef MPI
 636          call MPI_Finalize(ierr)
 637 #endif
 638          stop
 639         else if (mark.eq.0) then
 640           write (iout,*)&
 641           "ERROR: the inertia matrix is singular for chain",ichain
 642 #ifdef MPI
 643           call MPI_Finalize(ierr)
 644 #endif
 645         else if (mark.eq.1) then
 646           if (gmatout) then
 647           write (iout,*) "The transformed five-diagonal inertia matrix"
 648           write (iout,'(a,i5)') 'Chain',ichain
 649           do i=innt,inct
 650             if (i.lt.inct-1) then
 651               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
 652             else if (i.eq.inct-1) then
 653               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
 654             else
 655               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
 656             endif
 657           enddo
 658           endif
 659         endif
 660       enddo
 661 ! Diagonalization of the pentadiagonal matrix
 662 #ifdef TIMING
 663       time00=MPI_Wtime()
 664 #endif
 665
 666
 667 #endif
 668 !CRYSTBOND
 669 #else
 670       dimen=(nct-nnt+1)+nside
 671       dimen1=(nct-nnt)+(nct-nnt+1)
 672       dimen3=dimen*3
 673       write (iout,*) "nnt",nnt," nct",nct," nside",nside
 674 ! Old Gmatrix
 675 #ifdef MPI
 676       if (nfgtasks.gt.1) then
 677       time00=MPI_Wtime()
 678       call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
 679       time_Bcast=time_Bcast+MPI_Wtime()-time00
 680       call int_bounds(dimen,igmult_start,igmult_end)
 681       igmult_start=igmult_start-1
 682       call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,&
 683           ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 684       my_ng_count=igmult_end-igmult_start
 685       call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,&
 686           MPI_INTEGER,FG_COMM,IERROR)
 687       write (iout,*) 'Processor:',fg_rank,' CG group',kolor,&
 688        ' absolute rank',myrank,' igmult_start',igmult_start,&
 689        ' igmult_end',igmult_end,' count',my_ng_count
 690       write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 691       write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 692       call flush(iout)
 693       else
 694 #endif
 695       igmult_start=1
 696       igmult_end=dimen
 697       my_ng_count=dimen
 698 #ifdef MPI
 699       endif
 700 #endif
 701 !      write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
 702 !  Zeroing out A and fricmat
 703       do i=1,dimen
 704         do j=1,dimen
 705           A(i,j)=0.0D0
 706         enddo
 707       enddo
 708 !  Diagonal elements of the dC part of A and the respective friction coefficients
 709       ind=1
 710       ind1=0
 711 !      print *,"TUTUTUT?!",nnt,nct-1
 712       do i=nnt,nct-1
 713         mnum=molnum(i)
 714         ind=ind+1
 715         ind1=ind1+1
 716         coeff=0.25d0*IP(mnum)
 717         print *,"TU",i,itype(i,mnum),mnum
 718         if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
 719         print *,"TU2",i,coeff,mp(mnum)
 720         massvec(ind1)=mp(mnum)
 721         Gmat(ind,ind)=coeff
 722 !        print *,"i",mp(mnum)
 723         A(ind1,ind)=0.5d0
 724       enddo
 725
 726 !  Off-diagonal elements of the dC part of A
 727       k=3
 728       do i=1,nct-nnt
 729         do j=1,i
 730           A(i,j)=1.0d0
 731         enddo
 732       enddo
 733 !  Diagonal elements of the dX part of A and the respective friction coefficients
 734       m=nct-nnt
 735       m1=nct-nnt+1
 736       ind=0
 737       ind1=0
 738       do i=1,2
 739       msc(ntyp1_molec(i),i)=1.0d0
 740       enddo
 741       msc(ntyp1_molec(4),4)=1.0d0
 742 !      print *,"TU3",ntyp1_molec(4)-1
 743        msc(ntyp1_molec(4)-1,4)=1.0d0
 744       do i=nnt,nct
 745         ind=ind+1
 746         ii = ind+m
 747         mnum=molnum(i)
 748         iti=itype(i,mnum)
 749 !        if (mnum.eq.5) then
 750 !        mscab=0.0
 751 !        else
 752         mscab=msc(iabs(iti),mnum)
 753 !        endif
 754         massvec(ii)=mscab
 755
 756         if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.lt.4) then
 757           ind1=ind1+1
 758           ii1= ind1+m1
 759           A(ii,ii1)=1.0d0
 760           Gmat(ii1,ii1)=ISC(iabs(iti),mnum)
 761         endif
 762       enddo
 763 !  Off-diagonal elements of the dX part of A
 764       ind=0
 765       k=nct-nnt
 766       do i=nnt,nct
 767         mnum=molnum(i)
 768         iti=itype(i,mnum)
 769         ind=ind+1
 770         do j=nnt,i
 771           ii = ind
 772           jj = j-nnt+1
 773           A(k+ii,jj)=1.0d0
 774         enddo
 775       enddo
 776       if (lprn) then
 777         write (iout,*)
 778         write (iout,*) "Vector massvec"
 779         do i=1,dimen1
 780           write (iout,*) i,massvec(i)
 781         enddo
 782         write (iout,'(//a)') "A"
 783         call matout(dimen,dimen1,nres2,nres2,A)
 784       endif
 785
 786 ! Calculate the G matrix (store in Gmat)
 787       do k=1,dimen
 788        do i=1,dimen
 789          dtdi=0.0d0
 790          do j=1,dimen1
 791            dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
 792          enddo
 793          Gmat(k,i)=Gmat(k,i)+dtdi
 794        enddo
 795       enddo
 796
 797       if (lprn) then
 798         write (iout,'(//a)') "Gmat"
 799         call matout(dimen,dimen,nres2,nres2,Gmat)
 800       endif
 801       do i=1,dimen
 802         do j=1,dimen
 803           Ginv(i,j)=0.0d0
 804           Gcopy(i,j)=Gmat(i,j)
 805         enddo
 806         Ginv(i,i)=1.0d0
 807       enddo
 808 ! Invert the G matrix
 809       call MATINVERT(dimen,nres2,Gcopy,Ginv,osob)
 810       if (lprn) then
 811         write (iout,'(//a)') "Ginv"
 812         call matout(dimen,dimen,nres2,nres2,Ginv)
 813       endif
 814
 815       Bmat=0.0d0
 816       Bmat(1,1)=1.0d0
 817       j=2
 818       do i=nnt,nct
 819         mnum=molnum(i)
 820         if ((itype(i,1).eq.10).or.(itype(i,mnum).eq.ntyp1_molec(mnum))&
 821           .or.(mnum.ge.4)) then
 822           Bmat(i-nnt+2,j-1)=-1
 823           Bmat(i-nnt+2,j)=1
 824           j=j+1
 825         else
 826           if (i.lt.nct) then
 827             Bmat(i-nnt+2,j-1)=-1
 828             Bmat(i-nnt+2,j+1)=1
 829             Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
 830             Bmat(i-nnt+1+(nct-nnt)+1,j)=1
 831             j=j+2
 832           else
 833             Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
 834             Bmat(i-nnt+1+(nct-nnt)+1,j)=1
 835           endif
 836         endif
 837       enddo
 838       write (iout,*) "j",j," dimen",dimen
 839       if (lprn) then
 840         write (iout,'(//a)') "Bmat"
 841         call matout(dimen,dimen,nres2,nres2,Bmat)
 842       endif
 843       Gcopy=0.0d0
 844       write(iout,*) "before Gcopy",dimen,nres*2
 845 #ifdef TESTFIVEDIAG
 846       do i=1,dimen
 847         do j=1,dimen
 848 !            write (iout,*) "myij",i,j
 849             do k=1,dimen
 850                do l=1,dimen
 851 !                 write(iout,*) "i,j,k,l",i,j,k,l
 852                  Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j)
 853                enddo
 854             enddo
 855         enddo
 856       enddo
 857 #endif
 858       if (lprn) then
 859         write (iout,'(//a)') "Gmat-transformed"
 860         call matout(dimen,dimen,nres2,nres2,Gcopy)
 861       endif
 862 #ifdef MPI
 863       if (nfgtasks.gt.1) then
 864         myginv_ng_count=nres2*my_ng_count
 865         call MPI_Allgather(nres2*igmult_start,1,MPI_INTEGER,&
 866           nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 867         call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,&
 868           nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 869         write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 870         write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 871         call flush(iout)
 872 !        call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 873 !     &    nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
 874 !     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 875 !        call MPI_Barrier(FG_COMM,IERR)
 876         time00=MPI_Wtime()
 877         call MPI_Scatterv(ginv(1,1),nginv_counts(0),&
 878           nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),&
 879           myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 880 #ifdef TIMING
 881         time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 882 #endif
 883         do i=1,dimen
 884           do j=1,2*my_ng_count
 885             ginv(j,i)=gcopy(i,j)
 886           enddo
 887         enddo
 888 !        write (iout,*) "Master's chunk of ginv"
 889 !        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
 890       endif
 891 #endif
 892       if (osob) then
 893         write (iout,*) "The G matrix is singular."
 894         write (iout,'(//a)') "Gmat-transformed"
 895         call matout(dimen,dimen,nres2,nres2,Gcopy)
 896         stop
 897       endif
 898 ! Compute G**(-1/2) and G**(1/2)
 899       ind=0
 900       do i=1,dimen
 901         do j=1,i
 902           ind=ind+1
 903           Ghalf(ind)=Gmat(i,j)
 904         enddo
 905       enddo
 906       call gldiag(nres2,dimen,dimen,Ghalf,work,Geigen,Gvec,&
 907         ierr,iwork)
 908       if (lprn) then
 909         write (iout,'(//a)') &
 910          "Eigenvectors and eigenvalues of the G matrix"
 911         call eigout(dimen,dimen,nres2,nres2,Gvec,Geigen)
 912       endif
 913       do i=1,dimen
 914         sqreig(i)=dsqrt(Geigen(i))
 915       enddo
 916       do i=1,dimen
 917         do j=1,dimen
 918           Gsqrp(i,j)=0.0d0
 919           Gsqrm(i,j)=0.0d0
 920           Gcopy(i,j)=0.0d0
 921           do k=1,dimen
 922             Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
 923             Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
 924             Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
 925           enddo
 926         enddo
 927       enddo
 928       if (lprn) then
 929         write (iout,*) "Comparison of original and restored G"
 930         do i=1,dimen
 931           do j=1,dimen
 932             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),&
 933               Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
 934           enddo
 935         enddo
 936       endif
 937       if (allocated(Gcopy)) deallocate(Gcopy)
 938 #endif
 939 !write(iout,*) "end setup_MD_matr"
 940       return
 941       end subroutine setup_MD_matrices
 942 !-----------------------------------------------------------------------------
 943       subroutine EIGOUT(NC,NR,LM2,LM3,A,B)
 944
 945 !      implicit real*8 (a-h,o-z)
 946 !      include 'DIMENSIONS'
 947 !      include 'COMMON.IOUNITS'
 948       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 949       real(kind=8) :: A(LM2,LM3),B(LM2)
 950       KA=1
 951       KC=6
 952     1 KB=MIN0(KC,NC)
 953       WRITE(IOUT,600) (I,I=KA,KB)
 954       WRITE(IOUT,601) (B(I),I=KA,KB)
 955       WRITE(IOUT,602)
 956     2 N=0
 957       DO 3  I=1,NR
 958       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 959       N=N+1
 960       IF(N.LT.10) GO TO 3
 961       WRITE(IOUT,602)
 962       N=0
 963     3 CONTINUE
 964     4 IF (KB.EQ.NC) RETURN
 965       KA=KC+1
 966       KC=KC+6
 967       GO TO 1
 968   600 FORMAT (// 9H ROOT NO.,I4,9I11)
 969   601 FORMAT (/5X,10(1PE11.4))
 970   602 FORMAT (2H  )
 971   603 FORMAT (I5,10F11.5)
 972   604 FORMAT (1H1)
 973       end subroutine EIGOUT
 974 !-----------------------------------------------------------------------------
 975       subroutine MATOUT(NC,NR,LM2,LM3,A)
 976
 977 !      implicit real*8 (a-h,o-z)
 978 !      include 'DIMENSIONS'
 979 !      include 'COMMON.IOUNITS'
 980       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 981       real(kind=8) :: A(LM2,LM3)
 982       KA=1
 983       KC=6
 984     1 KB=MIN0(KC,NC)
 985       WRITE(IOUT,600) (I,I=KA,KB)
 986       WRITE(IOUT,602)
 987     2 N=0
 988       DO 3  I=1,NR
 989       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 990       N=N+1
 991       IF(N.LT.10) GO TO 3
 992       WRITE(IOUT,602)
 993       N=0
 994     3 CONTINUE
 995     4 IF (KB.EQ.NC) RETURN
 996       KA=KC+1
 997       KC=KC+6
 998       GO TO 1
 999   600 FORMAT (//5x,9I11)
1000   602 FORMAT (2H  )
1001   603 FORMAT (I5,10F11.3)
1002   604 FORMAT (1H1)
1003       end subroutine MATOUT
1004 !-----------------------------------------------------------------------------
1005       subroutine MATOUT1(NC,NR,LM2,LM3,A)
1006
1007 !      implicit real*8 (a-h,o-z)
1008 !      include 'DIMENSIONS'
1009 !      include 'COMMON.IOUNITS'
1010       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
1011       real(kind=8) :: A(LM2,LM3)
1012       KA=1
1013       KC=21
1014     1 KB=MIN0(KC,NC)
1015       WRITE(IOUT,600) (I,I=KA,KB)
1016       WRITE(IOUT,602)
1017     2 N=0
1018       DO 3  I=1,NR
1019       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
1020       N=N+1
1021       IF(N.LT.3) GO TO 3
1022       WRITE(IOUT,602)
1023       N=0
1024     3 CONTINUE
1025     4 IF (KB.EQ.NC) RETURN
1026       KA=KC+1
1027       KC=KC+21
1028       GO TO 1
1029   600 FORMAT (//5x,7(3I5,2x))
1030   602 FORMAT (2H  )
1031   603 FORMAT (I5,7(3F5.1,2x))
1032   604 FORMAT (1H1)
1033       end subroutine MATOUT1
1034 !-----------------------------------------------------------------------------
1035       subroutine MATOUT2(NC,NR,LM2,LM3,A)
1036
1037 !      implicit real*8 (a-h,o-z)
1038 !      include 'DIMENSIONS'
1039 !      include 'COMMON.IOUNITS'
1040       integer :: I,J,KA,KC,KB,N
1041       integer :: LM2,LM3,NC,NR
1042       real(kind=8) :: A(LM2,LM3)
1043       KA=1
1044       KC=12
1045     1 KB=MIN0(KC,NC)
1046       WRITE(IOUT,600) (I,I=KA,KB)
1047       WRITE(IOUT,602)
1048     2 N=0
1049       DO 3  I=1,NR
1050       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
1051       N=N+1
1052       IF(N.LT.3) GO TO 3
1053       WRITE(IOUT,602)
1054       N=0
1055     3 CONTINUE
1056     4 IF (KB.EQ.NC) RETURN
1057       KA=KC+1
1058       KC=KC+12
1059       GO TO 1
1060   600 FORMAT (//5x,4(3I9,2x))
1061   602 FORMAT (2H  )
1062   603 FORMAT (I5,4(3F9.3,2x))
1063   604 FORMAT (1H1)
1064       end subroutine MATOUT2
1065 #ifdef FIVEDIAG
1066       subroutine fivediagmult(n,DM,DU1,DU2,x,y)
1067       integer n
1068       double precision DM(n),DU1(n),DU2(n),x(n),y(n)
1069       integer i
1070       y(1)=DM(1)*x(1)+DU1(1)*x(2)+DU2(1)*x(3)
1071       y(2)=DU1(1)*x(1)+DM(2)*x(2)+DU1(2)*x(3)+DU2(2)*x(4)
1072       do i=3,n-2
1073         y(i)=DU2(i-2)*x(i-2)+DU1(i-1)*x(i-1)+DM(i)*x(i) &
1074            +DU1(i)*x(i+1)+DU2(i)*x(i+2)
1075       enddo
1076       y(n-1)=DU2(n-3)*x(n-3)+DU1(n-2)*x(n-2)+DM(n-1)*x(n-1) &
1077       +DU1(n-1)*x(n)
1078       y(n)=DU2(n-2)*x(n-2)+DU1(n-1)*x(n-1)+DM(n)*x(n)
1079       return
1080       end subroutine
1081 !c---------------------------------------------------------------------------
1082       subroutine fivediaginv_mult(ndim,forces,d_a_vec)
1083       use energy_data, only:nchain,chain_border,nct,nnt,molnum,&
1084       chain_border1,itype
1085       integer ndim
1086       double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim), &
1087        xsolv(ndim),d_a_vec(6*nres)
1088       integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev,mnum
1089       accel=0.0d0
1090       do j=1,3
1091 !Compute accelerations in Calpha and SC
1092         do ichain=1,nchain
1093           n=dimen_chain(ichain)
1094           iposc=iposd_chain(ichain)
1095           innt=chain_border(1,ichain)
1096           inct=chain_border(2,ichain)
1097           do i=iposc,iposc+n-1
1098             rs(i-iposc+1)=forces(3*(i-1)+j)
1099           enddo
1100 #ifdef DEBUG
1101           write (iout,*) "j",j," chain",ichain
1102           write (iout,*) "rs"
1103           write (iout,'(f10.5)') (rs(i),i=1,n)
1104 #endif
1105           call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
1106 #ifdef DEBUG
1107           write (iout,*) "xsolv"
1108           write (iout,'(f10.5)') (xsolv(i),i=1,n)
1109 #endif
1110           ind=1
1111           do i=innt,inct
1112             mnum=molnum(i)
1113             if (itype(i,1).eq.10.or.mnum.gt.2)then
1114               accel(j,i)=xsolv(ind)
1115               ind=ind+1
1116             else
1117               accel(j,i)=xsolv(ind)
1118               accel(j,i+nres)=xsolv(ind+1)
1119               ind=ind+2
1120             end if
1121           enddo
1122         enddo
1123       enddo
1124 !C Convert d_a to virtual-bon-vector basis
1125 #ifdef DEBUG
1126       write (iout,*) "accel in CA-SC basis"
1127       do i=1,nres
1128         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(accel(j,i),j=1,3),
1129      &      (accel(j,i+nres),j=1,3)
1130       enddo
1131       write (iout,*) "nnt",nnt
1132 #endif
1133       if (nnt.eq.1) then
1134         accel(:,0)=accel(:,1)
1135       endif
1136       do i=1,nres
1137       mnum=molnum(i)
1138         if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
1139         .or.mnum.ge.3) then
1140           do j=1,3
1141             accel(j,i)=accel(j,i+1)-accel(j,i)
1142           enddo
1143         else
1144           do j=1,3
1145             accel(j,i+nres)=accel(j,i+nres)-accel(j,i)
1146             accel(j,i)=accel(j,i+1)-accel(j,i)
1147           enddo
1148         end if
1149       enddo
1150       accel(:,nres)=0.0d0
1151       accel(:,nct)=0.0d0
1152       accel(:,2*nres)=0.0d0
1153       if (nnt.gt.1) then
1154         accel(:,0)=accel(:,1)
1155         accel(:,1)=0.0d0
1156       endif
1157       do ichain=2,nchain
1158         accel(:,chain_border1(1,ichain)-1)= &
1159          accel(:,chain_border1(1,ichain))
1160         accel(:,chain_border1(1,ichain))=0.0d0
1161       enddo
1162       do ichain=2,nchain
1163         accel(:,chain_border1(1,ichain)-1)= &
1164        accel(:,chain_border1(1,ichain)-1) &
1165         +accel(:,chain_border(2,ichain-1))
1166         accel(:,chain_border(2,ichain-1))=0.0d0
1167       enddo
1168 #ifdef DEBUG
1169       write (iout,*) "accel in fivediaginv_mult: 1"
1170       do i=0,2*nres
1171         write(iout,'(i5,3f10.5)') i,(accel(j,i),j=1,3)
1172       enddo
1173 #endif
1174       do j=1,3
1175         d_a_vec(j)=accel(j,0)
1176       enddo
1177       ind=3
1178       do i=nnt,nct-1
1179         do j=1,3
1180           d_a_vec(ind+j)=accel(j,i)
1181         enddo
1182         ind=ind+3
1183       enddo
1184       do i=nnt,nct
1185         mnum=molnum(i)
1186         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1187           .and.mnum.le.2) then
1188           do j=1,3
1189             d_a_vec(ind+j)=accel(j,i+nres)
1190           enddo
1191           ind=ind+3
1192         endif
1193       enddo
1194 #ifdef DEBUG
1195       write (iout,*) "d_a_vec"
1196       write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
1197 #endif
1198       return
1199       end subroutine
1200
1201
1202 #else
1203 !-----------------------------------------------------------------------------
1204       subroutine ginv_mult(z,d_a_tmp)
1205
1206       use geometry_data, only: nres
1207       use control_data
1208       use MPI_data
1209 !      implicit real*8 (a-h,o-z)
1210 !      include 'DIMENSIONS'
1211 #ifdef MPI
1212       include 'mpif.h'
1213       integer :: ierr,ierror
1214 #endif
1215 !      include 'COMMON.SETUP'
1216 !      include 'COMMON.TIME1'
1217 !      include 'COMMON.MD'
1218       real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
1219       real(kind=8),dimension(6*nres) :: temp    !(maxres6) maxres6=6*maxres
1220       real(kind=8) :: time00,time01
1221       integer :: i,j,k,ind
1222 #ifdef MPI
1223       if (nfgtasks.gt.1) then
1224         if (fg_rank.eq.0) then
1225 ! The matching BROADCAST for fg processors is called in ERGASTULUM
1226           time00=MPI_Wtime()
1227           call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
1228           time_Bcast=time_Bcast+MPI_Wtime()-time00
1229 !          print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
1230         endif
1231 !        write (2,*) "time00",time00
1232 !        write (2,*) "Before Scatterv"
1233 !        call flush(2)
1234 !        write (2,*) "Whole z (for FG master)"
1235 !        do i=1,dimen
1236 !          write (2,*) i,z(i)
1237 !        enddo
1238 !        call MPI_Barrier(FG_COMM,IERROR)
1239         time00=MPI_Wtime()
1240 !elwrite(iout,*) "do tej pory jest OK, MPI_Scatterv w ginv_mult"
1241         call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
1242           MPI_DOUBLE_PRECISION,&
1243           z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1244 !        write (2,*) "My chunk of z"
1245 !        do i=1,3*my_ng_count
1246 !          write (2,*) i,z(i)
1247 !        enddo
1248 !        write (2,*) "After SCATTERV"
1249 !        call flush(2)
1250 !        write (2,*) "MPI_Wtime",MPI_Wtime()
1251         time_scatter=time_scatter+MPI_Wtime()-time00
1252 #ifdef TIMING
1253         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
1254 #endif
1255 !        write (2,*) "time_scatter",time_scatter
1256 !        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
1257 !     &    my_ng_count
1258 !        call flush(2)
1259         time01=MPI_Wtime()
1260         do k=0,2
1261           do i=1,dimen
1262             ind=(i-1)*3+k+1
1263             temp(ind)=0.0d0
1264             do j=1,my_ng_count
1265 !              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
1266 !     &         Ginv(i,j),z((j-1)*3+k+1),
1267 !     &          Ginv(i,j)*z((j-1)*3+k+1)
1268 !              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
1269               temp(ind)=temp(ind)+Ginv(j,i)*z1((j-1)*3+k+1)
1270             enddo
1271           enddo
1272         enddo
1273         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1274 !        write (2,*) "Before REDUCE"
1275 !        call flush(2)
1276 !        write (2,*) "z before reduce"
1277 !        do i=1,dimen
1278 !          write (2,*) i,temp(i)
1279 !        enddo
1280         time00=MPI_Wtime()
1281         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
1282             MPI_SUM,king,FG_COMM,IERR)
1283         time_reduce=time_reduce+MPI_Wtime()-time00
1284 !        write (2,*) "After REDUCE"
1285 !        call flush(2)
1286       else
1287 #endif
1288 #ifdef TIMING
1289         time01=MPI_Wtime()
1290 #endif
1291         do k=0,2
1292           do i=1,dimen
1293             ind=(i-1)*3+k+1
1294             d_a_tmp(ind)=0.0d0
1295             do j=1,dimen
1296 !              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
1297 !              call flush(2)
1298 !     &         Ginv(i,j),z((j-1)*3+k+1),
1299 !     &          Ginv(i,j)*z((j-1)*3+k+1)
1300               d_a_tmp(ind)=d_a_tmp(ind) &
1301                                +Ginv(j,i)*z((j-1)*3+k+1)
1302 !              d_a_tmp(ind)=d_a_tmp(ind)
1303 !     &                         +Ginv(i,j)*z((j-1)*3+k+1)
1304             enddo
1305           enddo
1306         enddo
1307 #ifdef TIMING
1308         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1309 #endif
1310 #ifdef MPI
1311       endif
1312 #endif
1313       return
1314       end subroutine ginv_mult
1315 !-----------------------------------------------------------------------------
1316 #ifdef GINV_MULT
1317       subroutine ginv_mult_test(z,d_a_tmp)
1318
1319 !      include 'DIMENSIONS'
1320 !el      integer :: dimen
1321 !      include 'COMMON.MD'
1322       real(kind=8),dimension(dimen) :: z,d_a_tmp
1323       real(kind=8),dimension(dimen/3) :: ztmp,dtmp
1324       integer :: i,j,k,ind
1325 !      do i=1,dimen
1326 !        d_a_tmp(i)=0.0d0
1327 !        do j=1,dimen
1328 !          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
1329 !        enddo
1330 !      enddo
1331 !
1332 !      return
1333
1334 !ibm* unroll(3)
1335       do k=0,2
1336        do j=1,dimen/3
1337         ztmp(j)=z((j-1)*3+k+1)
1338        enddo
1339
1340        call alignx(16,ztmp(1))
1341        call alignx(16,dtmp(1))
1342        call alignx(16,Ginv(1,1))
1343
1344        do i=1,dimen/3
1345         dtmp(i)=0.0d0
1346         do j=1,dimen/3
1347            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
1348         enddo
1349        enddo
1350        do i=1,dimen/3
1351         ind=(i-1)*3+k+1
1352         d_a_tmp(ind)=dtmp(i)
1353        enddo
1354       enddo
1355       return
1356       end subroutine ginv_mult_test
1357 #endif
1358 !-----------------------------------------------------------------------------
1359       subroutine fricmat_mult(z,d_a_tmp)
1360
1361       use geometry_data, only: nres
1362       use control_data
1363       use MPI_data
1364 !      include 'DIMENSIONS'
1365 #ifdef MPI
1366       include 'mpif.h'
1367       integer :: IERROR,ierr
1368 #endif
1369 !      include 'COMMON.MD'
1370 !      include 'COMMON.IOUNITS'
1371 !      include 'COMMON.SETUP'
1372 !      include 'COMMON.TIME1'
1373 !#ifndef LANG0
1374 !      include 'COMMON.LANGEVIN'
1375 !#else
1376 !      include 'COMMON.LANGEVIN.lang0'
1377 !#endif
1378       real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
1379       real(kind=8),dimension(6*nres) :: temp    !(maxres6) maxres6=6*maxres
1380       real(kind=8) :: time00,time01
1381       integer :: i,j,k,ind,nres2
1382       nres2=2*nres
1383 !el      if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
1384
1385 #ifdef MPI
1386       if (nfgtasks.gt.1) then
1387         if (fg_rank.eq.0) then
1388 ! The matching BROADCAST for fg processors is called in ERGASTULUM
1389           time00=MPI_Wtime()
1390           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
1391           time_Bcast=time_Bcast+MPI_Wtime()-time00
1392 !          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
1393         endif
1394 !        call MPI_Barrier(FG_COMM,IERROR)
1395         time00=MPI_Wtime()
1396         call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
1397           MPI_DOUBLE_PRECISION,&
1398           z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1399 !        write (2,*) "My chunk of z"
1400 !        do i=1,3*my_ng_count
1401 !          write (2,*) i,z(i)
1402 !        enddo
1403         time_scatter=time_scatter+MPI_Wtime()-time00
1404 #ifdef TIMING
1405         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
1406 #endif
1407         time01=MPI_Wtime()
1408         do k=0,2
1409           do i=1,dimen
1410             ind=(i-1)*3+k+1
1411             temp(ind)=0.0d0
1412             do j=1,my_ng_count
1413               temp(ind)=temp(ind)-fricmat(j,i)*z1((j-1)*3+k+1)
1414             enddo
1415           enddo
1416         enddo
1417         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1418 !        write (2,*) "Before REDUCE"
1419 !        write (2,*) "d_a_tmp before reduce"
1420 !        do i=1,dimen3
1421 !          write (2,*) i,temp(i)
1422 !        enddo
1423 !        call flush(2)
1424         time00=MPI_Wtime()
1425         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
1426             MPI_SUM,king,FG_COMM,IERR)
1427         time_reduce=time_reduce+MPI_Wtime()-time00
1428 !        write (2,*) "After REDUCE"
1429 !        call flush(2)
1430       else
1431 #endif
1432 #ifdef TIMING
1433         time01=MPI_Wtime()
1434 #endif
1435         do k=0,2
1436          do i=1,dimen
1437           ind=(i-1)*3+k+1
1438           d_a_tmp(ind)=0.0d0
1439           do j=1,dimen
1440              d_a_tmp(ind)=d_a_tmp(ind) &
1441                                  -fricmat(j,i)*z((j-1)*3+k+1)
1442           enddo
1443          enddo
1444         enddo
1445 #ifdef TIMING
1446         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1447 #endif
1448 #ifdef MPI
1449       endif
1450 #endif
1451 !      write (iout,*) "Vector d_a"
1452 !      do i=1,dimen3
1453 !        write (2,*) i,d_a_tmp(i)
1454 !      enddo
1455       return
1456       end subroutine fricmat_mult
1457 !-----------------------------------------------------------------------------
1458 #endif
1459 !-----------------------------------------------------------------------------
1460       end module REMD