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       use geometry_data, only: nres
  28       use control_data    !el, only: mucadyn,lmuca
  29 #ifdef MPI
  30        include 'mpif.h'
  31       real(kind=8) :: time00
  32 #endif
  33 !       include 'COMMON.VAR'
  34 !       include 'COMMON.CHAIN'
  35 !       include 'COMMON.DERIV'
  36 !       include 'COMMON.GEO'
  37 !       include 'COMMON.LOCAL'
  38 !       include 'COMMON.INTERACT'
  39 !       include 'COMMON.MD'
  40 !       include 'COMMON.IOUNITS'
  41 !       include 'COMMON.CONTROL'
  42 !       include 'COMMON.MUCA'
  43 !      include 'COMMON.TIME1'
  44        integer ::i,j,ind,itime,mnum
  45        real(kind=8) :: zapas(6*nres) !,muca_factor      !maxres6=6*maxres
  46        real(kind=8) :: rs(dimen),xsolv(dimen)
  47 #ifdef CHECK5DSOL
  48        real(kind=8) :: rscheck(dimen),rsold(dimen)
  49 #endif
  50        logical :: lprn = .true.
  51 !el       common /cipiszcze/ itime
  52        itime = itt_comm
  53
  54 #ifdef TIMING
  55        time00=MPI_Wtime()
  56 #endif
  57 #ifdef FIVEDIAG
  58       if (lprn) then
  59         write (iout,*) "Potential forces backbone"
  60         do i=nnt,nct
  61           write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3)
  62         enddo
  63         write (iout,*) "Potential forces sidechain"
  64         do i=nnt,nct
  65           if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) &
  66             write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
  67         enddo
  68       endif
  69        do j=1,3
  70          ind=1
  71          do i=nnt,nct
  72            if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)then
  73             rs(ind)=-gcart(j,i)-gxcart(j,i)
  74             ind=ind+1
  75            else
  76             rs(ind)=-gcart(j,i)
  77             rs(ind+1)=-gxcart(j,i)
  78             ind=ind+2
  79            end if
  80          enddo
  81 #ifdef CHECK5DSOL
  82          rsold=rs
  83 #endif
  84          if (lprn) then
  85            write(iout,*) "RHS of the 5-diag equations system"
  86            do i=1,dimen
  87              write(iout,*) i,rs(i)
  88            enddo
  89          endif
  90
  91          call FDISYS (dimen,DM,DU1,DU2,rs,xsolv)
  92          if (lprn) then
  93            write (iout,*) "Solution of the 5-diagonal equations system"
  94            do i=1,dimen
  95              write (iout,'(i5,f10.5)') i,xsolv(i)
  96            enddo
  97          endif
  98 #ifdef CHECK5DSOL
  99 ! Check the solution
 100          rscheck(1)=DMorig(1)*xsolv(1)+DU1orig(1)*xsolv(2)+&
 101            DU2orig(1)*xsolv(3)
 102          rscheck(2)=DU1orig(1)*xsolv(1)+DMorig(2)*xsolv(2)+&
 103            DU1orig(2)*xsolv(3)+DU2orig(2)*xsolv(4)
 104
 105          do i=3,dimen-2
 106           rscheck(i)=DU2orig(i-2)*xsolv(i-2)+DU1orig(i-1)*&
 107             xsolv(i-1)+DMorig(i)*xsolv(i)+DU1orig(i)*&
 108              xsolv(i+1)+DU2orig(i)*xsolv(i+2)
 109          enddo
 110        rscheck(dimen-1)=DU2orig(dimen-3)*xsolv(dimen-3)+&
 111          DU1orig(dimen-2)*xsolv(dimen-2)+DMorig(dimen-1)*&
 112           xsolv(dimen-1)+DU1orig(dimen-1)*&
 113            xsolv(dimen)
 114        rscheck(dimen)=DU2orig(dimen-2)*xsolv(dimen-2)+DU1orig(dimen-1)*&
 115           xsolv(dimen-1)+DMorig(dimen)*xsolv(dimen)
 116
 117          do i=1,dimen
 118             write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),&
 119             "ratio",rscheck(i)/rsold(i)
 120            enddo
 121 ! end check
 122 #endif
 123          ind=1
 124          do i=nnt,nct
 125           mnum=molnum(i)
 126           if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum))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
 136        if (lprn) then
 137        do i=nnt,nct
 138         write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
 139        enddo
 140        do i=nnt,nct
 141         write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
 142        enddo
 143        endif
 144        do j=1,3
 145          d_a(j,0)=d_a(j,nnt)
 146        enddo
 147        do i=nnt,nct
 148           mnum=molnum(i)
 149 !         if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1) then
 150           if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum))then
 151            do j=1,3
 152              d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 153            enddo
 154          else
 155            do j=1,3
 156              d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
 157              d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 158            enddo
 159
 160
 161          end if
 162        enddo
 163
 164       if (lprn) then
 165         write(iout,*) 'acceleration 3D FIVEDIAG'
 166         write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
 167         do i=nnt,nct-1
 168          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 169         enddo
 170         do i=nnt,nct
 171          write (iout,'(i3,3f10.5,3x,3f10.5)') &
 172            i+nres,(d_a(j,i+nres),j=1,3)
 173         enddo
 174       endif
 175 #else
 176 ! Old procedure
 177        do j=1,3
 178          zapas(j)=-gcart(j,0)
 179        enddo
 180       ind=3
 181       if (lprn) then
 182         write (iout,*) "Potential forces backbone"
 183       endif
 184       do i=nnt,nct-1
 185         if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
 186           i,(-gcart(j,i),j=1,3)
 187         do j=1,3
 188           ind=ind+1
 189           zapas(ind)=-gcart(j,i)
 190         enddo
 191       enddo
 192       if (lprn) write (iout,*) "Potential forces sidechain"
 193       do i=nnt,nct
 194           mnum=molnum(i)
 195         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
 196           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
 197              i,(-gxcart(j,i),j=1,3)
 198           do j=1,3
 199             ind=ind+1
 200             zapas(ind)=-gxcart(j,i)
 201           enddo
 202         endif
 203       enddo
 204
 205       call ginv_mult(zapas,d_a_work)
 206
 207       do j=1,3
 208         d_a(j,0)=d_a_work(j)
 209       enddo
 210       ind=3
 211       do i=nnt,nct-1
 212         do j=1,3
 213           ind=ind+1
 214           d_a(j,i)=d_a_work(ind)
 215         enddo
 216       enddo
 217       do i=nnt,nct
 218        mnum=molnum(i)
 219 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
 220          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
 221           do j=1,3
 222             ind=ind+1
 223             d_a(j,i+nres)=d_a_work(ind)
 224           enddo
 225         endif
 226       enddo
 227
 228 #endif
 229       if(lmuca) then
 230        imtime=imtime+1
 231        if(mucadyn.gt.0) call muca_update(potE)
 232        factor=muca_factor(potE)*t_bath*Rb
 233
 234 !d       print *,'lmuca ',factor,potE
 235        do j=1,3
 236           d_a(j,0)=d_a(j,0)*factor
 237        enddo
 238        do i=nnt,nct-1
 239          do j=1,3
 240           d_a(j,i)=d_a(j,i)*factor
 241          enddo
 242        enddo
 243        do i=nnt,nct
 244          do j=1,3
 245           d_a(j,i+nres)=d_a(j,i+nres)*factor
 246          enddo
 247        enddo
 248
 249       endif
 250
 251       if (lprn) then
 252         write(iout,*) 'acceleration 3D'
 253         write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
 254         do i=nnt,nct-1
 255          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 256         enddo
 257         do i=nnt,nct
 258          write (iout,'(i3,3f10.5,3x,3f10.5)') &
 259            i+nres,(d_a(j,i+nres),j=1,3)
 260         enddo
 261       endif
 262 #ifdef TIMING
 263       time_lagrangian=time_lagrangian+MPI_Wtime()-time00
 264 #endif
 265       return
 266       end subroutine lagrangian
 267 !-----------------------------------------------------------------------------
 268       subroutine setup_MD_matrices
 269
 270       use geometry_data, only: nres,nside
 271       use control_data
 272       use MPI_data
 273       use energy_data
 274       use geometry, only:int_bounds
 275       use md_calc
 276 !      implicit real*8 (a-h,o-z)
 277 !      include 'DIMENSIONS'
 278 #ifdef MPI
 279       include 'mpif.h'
 280       integer :: ierror
 281       real(kind=8) :: time00
 282       real(kind=8) ,allocatable, dimension(:)  :: DDM,DDU1,DDU2
 283 #endif
 284 !      include 'COMMON.SETUP'
 285 !      include 'COMMON.VAR'
 286 !      include 'COMMON.CHAIN'
 287 !      include 'COMMON.DERIV'
 288 !      include 'COMMON.GEO'
 289 !      include 'COMMON.LOCAL'
 290 !      include 'COMMON.INTERACT'
 291 !      include 'COMMON.MD'
 292 !#ifndef LANG0
 293 !      include 'COMMON.LANGEVIN'
 294 !#else
 295 !      include 'COMMON.LANGEVIN.lang0'
 296 !#endif
 297 !      include 'COMMON.IOUNITS'
 298 !      include 'COMMON.TIME1'
 299       logical :: lprn = .true.
 300       logical :: osob
 301 #ifndef FIVEDIAG
 302       real(kind=8),dimension(2*nres,2*nres) :: Bmat,matmult
 303 #endif
 304       real(kind=8) :: dtdi
 305       real(kind=8),dimension(2*nres) :: massvec,sqreig  !(maxres2) maxres2=2*maxres
 306       real(kind=8) :: relfeh,eps1,eps2
 307 !el      real(kind=8),dimension(:),allocatable :: Ghalf
 308 !el      real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf   !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
 309 !el      real(kind=8),dimension(2*nres,2*nres) :: Gcopy !(maxres2,maxres2)
 310 !el      real(kind=8),dimension(:,:),allocatable :: Gcopy
 311       real(kind=8),dimension(8*6*nres) :: work  !(8*maxres6)
 312       integer,dimension(6*nres) :: iwork        !(maxres6) maxres6=6*maxres
 313 !el      common /przechowalnia/ Gcopy,Ghalf
 314       real(kind=8) :: coeff
 315       integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
 316       real(kind=8) :: ip4
 317       integer :: iz,mnum
 318       relfeh=1.0d-14
 319       nres2=2*nres
 320
 321       write (iout,*) "before FIVEDIAG"
 322 #ifndef FIVEDIAG
 323       write (iout,*) "ALLOCATE"
 324       if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
 325       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
 326 #endif
 327 !
 328 ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
 329 ! inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
 330 !
 331 ! Determine the number of degrees of freedom (dimen) and the number of
 332 ! sites (dimen1)
 333       dimen=(nct-nnt+1)+nside
 334       dimen1=(nct-nnt)+(nct-nnt+1)
 335       dimen3=dimen*3
 336       write (iout,*) "nnt",nnt," nct",nct," nside",nside
 337 #ifdef FIVEDIAG
 338        ip4=ip/4
 339        DM(1)=mp/4+ip4
 340        if (iabs(itype(nnt)).eq.10) then
 341          DM(1)=DM(1)+msc(10)
 342          ind=2
 343          nind=1
 344        else
 345          DM(1)=DM(1)+isc(iabs(itype(nnt)))
 346          DM(2)=msc(iabs(itype(nnt)))+isc(iabs(itype(nnt)))
 347          ind=3
 348          nind=2
 349        endif
 350        do i=nnt+1,nct-1
 351 !         if (iabs(itype(i,1)).eq.ntyp1) cycle
 352          DM(ind)=2*ip4+mp/2
 353          if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) then
 354            if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10)
 355            ind=ind+1
 356          else
 357            DM(ind)=DM(ind)+isc(iabs(itype(i,1)))
 358            DM(ind+1)=msc(iabs(itype(i,1)))+isc(iabs(itype(i,1)))
 359            ind=ind+2
 360          endif
 361        enddo
 362        if (nct.gt.nnt) then
 363        DM(ind)=ip4+mp/4
 364        if (iabs(itype(nct)).eq.10) then
 365          DM(ind)=DM(ind)+msc(10)
 366          nind=ind
 367        else
 368          DM(ind)=DM(ind)+isc(iabs(itype(nct)))
 369          DM(ind+1)=msc(iabs(itype(nct)))+isc(iabs(itype(nct)))
 370          nind=ind+1
 371         endif
 372         endif
 373
 374
 375         ind=1
 376         do i=nnt,nct
 377         if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
 378           DU1(ind)=-isc(iabs(itype(i,1)))
 379           DU1(ind+1)=0.0d0
 380           ind=ind+2
 381         else
 382           DU1(ind)=mp/4-ip4
 383           ind=ind+1
 384         endif
 385         enddo
 386
 387        ind=1
 388        do i=nnt,nct-1
 389        mnum=molnum(i)
 390 !        if (iabs(itype(i,1)).eq.ntyp1) cycle
 391         write (iout,*) "i",i," itype",itype(i,1),ntyp1
 392         if (iabs(itype(i,1)).ne.10 .and. &
 393           iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
 394          DU2(ind)=mp/4-ip4
 395          DU2(ind+1)=0.0d0
 396          ind=ind+2
 397         else
 398          DU2(ind)=0.0d0
 399          DU2(ind+1)=0.0d0
 400          ind=ind+1
 401         endif
 402        enddo
 403        DMorig=DM
 404        DU1orig=DU1
 405        DU2orig=DU2
 406        write (iout,*) "nind",nind," dimen",dimen
 407        if (nind.ne.dimen) then
 408          write (iout,*) "ERROR, dimensions differ"
 409 #ifdef MPI
 410          call MPI_Finalize(ierr)
 411 #endif
 412          stop
 413        endif
 414      write (iout,*) "The upper part of the five-diagonal inertia matrix"
 415        do i=1,dimen
 416          if (i.lt.dimen-1) then
 417            write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
 418          else if (i.eq.dimen-1) then
 419            write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
 420          else
 421            write (iout,'(i3,3f10.5)') i,DM(i)
 422          endif
 423        enddo
 424
 425         call FDISYP (dimen, DM, DU1, DU2, MARK)
 426
 427         if (mark.eq.-1) then
 428      write (iout,*) "ERROR: the inertia matrix is not positive definite"
 429 #ifdef MPI
 430          call MPI_Finalize(ierr)
 431 #endif
 432          stop
 433         else if (mark.eq.0) then
 434      write (iout,*) "ERROR: the inertia matrix is singular"
 435 #ifdef MPI
 436          call MPI_Finalize(ierr)
 437 #endif
 438         else if (mark.eq.1) then
 439           write (iout,*) "The transformed inertia matrix"
 440           do i=1,dimen
 441             if (i.lt.dimen-1) then
 442               write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
 443             else if (i.eq.dimen-1) then
 444               write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
 445             else
 446               write (iout,'(i3,3f10.5)') i,DM(i)
 447             endif
 448           enddo
 449         endif
 450 ! Diagonalization of the pentadiagonal matrix
 451         allocate(DDU1(2*nres))
 452         allocate(DDU2(2*nres))
 453         allocate(DDM(2*nres))
 454         do i=1,dimen-1
 455           DDU1(i+1)=DU1orig(i)
 456         enddo
 457         do i=1,dimen-2
 458           DDU2(i+2)=DU2orig(i)
 459         enddo
 460         DDM=DMorig
 461         call quindibisect2(DDM,DDU1,DDU2,dimen,1,dimen,eps1,relfeh,eps2,iz,geigen)
 462         if (lprn) then
 463           write (iout,*) &
 464            "Eigenvalues of the five-diagonal inertia matrix"
 465           do i=1,dimen
 466             write (iout,'(i5,f10.5)') i,geigen(i)
 467           enddo
 468         endif
 469         deallocate(DDU1)
 470         deallocate(DDU2)
 471         deallocate(DDM)
 472 #else
 473 ! Old Gmatrix
 474 #ifdef MPI
 475       if (nfgtasks.gt.1) then
 476       time00=MPI_Wtime()
 477       call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
 478       time_Bcast=time_Bcast+MPI_Wtime()-time00
 479       call int_bounds(dimen,igmult_start,igmult_end)
 480       igmult_start=igmult_start-1
 481       call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,&
 482           ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 483       my_ng_count=igmult_end-igmult_start
 484       call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,&
 485           MPI_INTEGER,FG_COMM,IERROR)
 486       write (iout,*) 'Processor:',fg_rank,' CG group',kolor,&
 487        ' absolute rank',myrank,' igmult_start',igmult_start,&
 488        ' igmult_end',igmult_end,' count',my_ng_count
 489       write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 490       write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 491       call flush(iout)
 492       else
 493 #endif
 494       igmult_start=1
 495       igmult_end=dimen
 496       my_ng_count=dimen
 497 #ifdef MPI
 498       endif
 499 #endif
 500 !      write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
 501 !  Zeroing out A and fricmat
 502       do i=1,dimen
 503         do j=1,dimen
 504           A(i,j)=0.0D0
 505         enddo
 506       enddo
 507 !  Diagonal elements of the dC part of A and the respective friction coefficients
 508       ind=1
 509       ind1=0
 510       print *,"TUTUTUT?!",nnt,nct-1
 511       do i=nnt,nct-1
 512         mnum=molnum(i)
 513         ind=ind+1
 514         ind1=ind1+1
 515         coeff=0.25d0*IP(mnum)
 516         massvec(ind1)=mp(mnum)
 517         Gmat(ind,ind)=coeff
 518         print *,"i",mp(mnum)
 519         A(ind1,ind)=0.5d0
 520       enddo
 521
 522 !  Off-diagonal elements of the dC part of A
 523       k=3
 524       do i=1,nct-nnt
 525         do j=1,i
 526           A(i,j)=1.0d0
 527         enddo
 528       enddo
 529 !  Diagonal elements of the dX part of A and the respective friction coefficients
 530       m=nct-nnt
 531       m1=nct-nnt+1
 532       ind=0
 533       ind1=0
 534       do i=1,2
 535       msc(ntyp1_molec(i),i)=1.0d0
 536       enddo
 537       do i=nnt,nct
 538         ind=ind+1
 539         ii = ind+m
 540         mnum=molnum(i)
 541         iti=itype(i,mnum)
 542         massvec(ii)=msc(iabs(iti),mnum)
 543         if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
 544           ind1=ind1+1
 545           ii1= ind1+m1
 546           A(ii,ii1)=1.0d0
 547           Gmat(ii1,ii1)=ISC(iabs(iti),mnum)
 548         endif
 549       enddo
 550 !  Off-diagonal elements of the dX part of A
 551       ind=0
 552       k=nct-nnt
 553       do i=nnt,nct
 554         iti=itype(i,1)
 555         ind=ind+1
 556         do j=nnt,i
 557           ii = ind
 558           jj = j-nnt+1
 559           A(k+ii,jj)=1.0d0
 560         enddo
 561       enddo
 562       if (lprn) then
 563         write (iout,*)
 564         write (iout,*) "Vector massvec"
 565         do i=1,dimen1
 566           write (iout,*) i,massvec(i)
 567         enddo
 568         write (iout,'(//a)') "A"
 569         call matout(dimen,dimen1,nres2,nres2,A)
 570       endif
 571
 572 ! Calculate the G matrix (store in Gmat)
 573       do k=1,dimen
 574        do i=1,dimen
 575          dtdi=0.0d0
 576          do j=1,dimen1
 577            dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
 578          enddo
 579          Gmat(k,i)=Gmat(k,i)+dtdi
 580        enddo
 581       enddo
 582
 583       if (lprn) then
 584         write (iout,'(//a)') "Gmat"
 585         call matout(dimen,dimen,nres2,nres2,Gmat)
 586       endif
 587       do i=1,dimen
 588         do j=1,dimen
 589           Ginv(i,j)=0.0d0
 590           Gcopy(i,j)=Gmat(i,j)
 591         enddo
 592         Ginv(i,i)=1.0d0
 593       enddo
 594 ! Invert the G matrix
 595       call MATINVERT(dimen,nres2,Gcopy,Ginv,osob)
 596       if (lprn) then
 597         write (iout,'(//a)') "Ginv"
 598         call matout(dimen,dimen,nres2,nres2,Ginv)
 599       endif
 600
 601       Bmat=0.0d0
 602       Bmat(1,1)=1.0d0
 603       j=2
 604       do i=nnt,nct
 605         mnum=molnum(i)
 606         if ((itype(i,1).eq.10).or.(itype(i,mnum).eq.ntyp1_molec(mnum))&
 607           .or.(mnum.eq.5)) then
 608           Bmat(i-nnt+2,j-1)=-1
 609           Bmat(i-nnt+2,j)=1
 610           j=j+1
 611         else
 612           if (i.lt.nct) then
 613             Bmat(i-nnt+2,j-1)=-1
 614             Bmat(i-nnt+2,j+1)=1
 615             Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
 616             Bmat(i-nnt+1+(nct-nnt)+1,j)=1
 617             j=j+2
 618           else
 619             Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
 620             Bmat(i-nnt+1+(nct-nnt)+1,j)=1
 621           endif
 622         endif
 623       enddo
 624       write (iout,*) "j",j," dimen",dimen
 625       if (lprn) then
 626         write (iout,'(//a)') "Bmat"
 627         call matout(dimen,dimen,nres2,nres2,Bmat)
 628       endif
 629       Gcopy=0.0d0
 630       do i=1,dimen
 631         do j=1,dimen
 632             do k=1,dimen
 633                do l=1,dimen
 634                  Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j)
 635                enddo
 636             enddo
 637         enddo
 638       enddo
 639       if (lprn) then
 640         write (iout,'(//a)') "Gmat-transformed"
 641         call matout(dimen,dimen,nres2,nres2,Gcopy)
 642       endif
 643 #ifdef MPI
 644       if (nfgtasks.gt.1) then
 645         myginv_ng_count=nres2*my_ng_count
 646         call MPI_Allgather(nres2*igmult_start,1,MPI_INTEGER,&
 647           nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 648         call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,&
 649           nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 650         write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 651         write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 652         call flush(iout)
 653 !        call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 654 !     &    nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
 655 !     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 656 !        call MPI_Barrier(FG_COMM,IERR)
 657         time00=MPI_Wtime()
 658         call MPI_Scatterv(ginv(1,1),nginv_counts(0),&
 659           nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),&
 660           myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 661 #ifdef TIMING
 662         time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 663 #endif
 664         do i=1,dimen
 665           do j=1,2*my_ng_count
 666             ginv(j,i)=gcopy(i,j)
 667           enddo
 668         enddo
 669 !        write (iout,*) "Master's chunk of ginv"
 670 !        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
 671       endif
 672 #endif
 673       if (osob) then
 674         write (iout,*) "The G matrix is singular."
 675         stop
 676       endif
 677 ! Compute G**(-1/2) and G**(1/2)
 678       ind=0
 679       do i=1,dimen
 680         do j=1,i
 681           ind=ind+1
 682           Ghalf(ind)=Gmat(i,j)
 683         enddo
 684       enddo
 685       call gldiag(nres2,dimen,dimen,Ghalf,work,Geigen,Gvec,&
 686         ierr,iwork)
 687       if (lprn) then
 688         write (iout,'(//a)') &
 689          "Eigenvectors and eigenvalues of the G matrix"
 690         call eigout(dimen,dimen,nres2,nres2,Gvec,Geigen)
 691       endif
 692       do i=1,dimen
 693         sqreig(i)=dsqrt(Geigen(i))
 694       enddo
 695       do i=1,dimen
 696         do j=1,dimen
 697           Gsqrp(i,j)=0.0d0
 698           Gsqrm(i,j)=0.0d0
 699           Gcopy(i,j)=0.0d0
 700           do k=1,dimen
 701             Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
 702             Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
 703             Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
 704           enddo
 705         enddo
 706       enddo
 707       if (lprn) then
 708         write (iout,*) "Comparison of original and restored G"
 709         do i=1,dimen
 710           do j=1,dimen
 711             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),&
 712               Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
 713           enddo
 714         enddo
 715       endif
 716       deallocate(Gcopy)
 717 #endif
 718       return
 719       end subroutine setup_MD_matrices
 720 !-----------------------------------------------------------------------------
 721       subroutine EIGOUT(NC,NR,LM2,LM3,A,B)
 722
 723 !      implicit real*8 (a-h,o-z)
 724 !      include 'DIMENSIONS'
 725 !      include 'COMMON.IOUNITS'
 726       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 727       real(kind=8) :: A(LM2,LM3),B(LM2)
 728       KA=1
 729       KC=6
 730     1 KB=MIN0(KC,NC)
 731       WRITE(IOUT,600) (I,I=KA,KB)
 732       WRITE(IOUT,601) (B(I),I=KA,KB)
 733       WRITE(IOUT,602)
 734     2 N=0
 735       DO 3  I=1,NR
 736       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 737       N=N+1
 738       IF(N.LT.10) GO TO 3
 739       WRITE(IOUT,602)
 740       N=0
 741     3 CONTINUE
 742     4 IF (KB.EQ.NC) RETURN
 743       KA=KC+1
 744       KC=KC+6
 745       GO TO 1
 746   600 FORMAT (// 9H ROOT NO.,I4,9I11)
 747   601 FORMAT (/5X,10(1PE11.4))
 748   602 FORMAT (2H  )
 749   603 FORMAT (I5,10F11.5)
 750   604 FORMAT (1H1)
 751       end subroutine EIGOUT
 752 !-----------------------------------------------------------------------------
 753       subroutine MATOUT(NC,NR,LM2,LM3,A)
 754
 755 !      implicit real*8 (a-h,o-z)
 756 !      include 'DIMENSIONS'
 757 !      include 'COMMON.IOUNITS'
 758       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 759       real(kind=8) :: A(LM2,LM3)
 760       KA=1
 761       KC=6
 762     1 KB=MIN0(KC,NC)
 763       WRITE(IOUT,600) (I,I=KA,KB)
 764       WRITE(IOUT,602)
 765     2 N=0
 766       DO 3  I=1,NR
 767       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 768       N=N+1
 769       IF(N.LT.10) GO TO 3
 770       WRITE(IOUT,602)
 771       N=0
 772     3 CONTINUE
 773     4 IF (KB.EQ.NC) RETURN
 774       KA=KC+1
 775       KC=KC+6
 776       GO TO 1
 777   600 FORMAT (//5x,9I11)
 778   602 FORMAT (2H  )
 779   603 FORMAT (I5,10F11.3)
 780   604 FORMAT (1H1)
 781       end subroutine MATOUT
 782 !-----------------------------------------------------------------------------
 783       subroutine MATOUT1(NC,NR,LM2,LM3,A)
 784
 785 !      implicit real*8 (a-h,o-z)
 786 !      include 'DIMENSIONS'
 787 !      include 'COMMON.IOUNITS'
 788       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 789       real(kind=8) :: A(LM2,LM3)
 790       KA=1
 791       KC=21
 792     1 KB=MIN0(KC,NC)
 793       WRITE(IOUT,600) (I,I=KA,KB)
 794       WRITE(IOUT,602)
 795     2 N=0
 796       DO 3  I=1,NR
 797       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 798       N=N+1
 799       IF(N.LT.3) GO TO 3
 800       WRITE(IOUT,602)
 801       N=0
 802     3 CONTINUE
 803     4 IF (KB.EQ.NC) RETURN
 804       KA=KC+1
 805       KC=KC+21
 806       GO TO 1
 807   600 FORMAT (//5x,7(3I5,2x))
 808   602 FORMAT (2H  )
 809   603 FORMAT (I5,7(3F5.1,2x))
 810   604 FORMAT (1H1)
 811       end subroutine MATOUT1
 812 !-----------------------------------------------------------------------------
 813       subroutine MATOUT2(NC,NR,LM2,LM3,A)
 814
 815 !      implicit real*8 (a-h,o-z)
 816 !      include 'DIMENSIONS'
 817 !      include 'COMMON.IOUNITS'
 818       integer :: I,J,KA,KC,KB,N
 819       integer :: LM2,LM3,NC,NR
 820       real(kind=8) :: A(LM2,LM3)
 821       KA=1
 822       KC=12
 823     1 KB=MIN0(KC,NC)
 824       WRITE(IOUT,600) (I,I=KA,KB)
 825       WRITE(IOUT,602)
 826     2 N=0
 827       DO 3  I=1,NR
 828       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 829       N=N+1
 830       IF(N.LT.3) GO TO 3
 831       WRITE(IOUT,602)
 832       N=0
 833     3 CONTINUE
 834     4 IF (KB.EQ.NC) RETURN
 835       KA=KC+1
 836       KC=KC+12
 837       GO TO 1
 838   600 FORMAT (//5x,4(3I9,2x))
 839   602 FORMAT (2H  )
 840   603 FORMAT (I5,4(3F9.3,2x))
 841   604 FORMAT (1H1)
 842       end subroutine MATOUT2
 843 !-----------------------------------------------------------------------------
 844       subroutine ginv_mult(z,d_a_tmp)
 845
 846       use geometry_data, only: nres
 847       use control_data
 848       use MPI_data
 849 !      implicit real*8 (a-h,o-z)
 850 !      include 'DIMENSIONS'
 851 #ifdef MPI
 852       include 'mpif.h'
 853       integer :: ierr,ierror
 854 #endif
 855 !      include 'COMMON.SETUP'
 856 !      include 'COMMON.TIME1'
 857 !      include 'COMMON.MD'
 858       real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
 859       real(kind=8),dimension(6*nres) :: temp    !(maxres6) maxres6=6*maxres
 860       real(kind=8) :: time00,time01
 861       integer :: i,j,k,ind
 862 #ifdef MPI
 863       if (nfgtasks.gt.1) then
 864         if (fg_rank.eq.0) then
 865 ! The matching BROADCAST for fg processors is called in ERGASTULUM
 866           time00=MPI_Wtime()
 867           call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
 868           time_Bcast=time_Bcast+MPI_Wtime()-time00
 869 !          print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
 870         endif
 871 !        write (2,*) "time00",time00
 872 !        write (2,*) "Before Scatterv"
 873 !        call flush(2)
 874 !        write (2,*) "Whole z (for FG master)"
 875 !        do i=1,dimen
 876 !          write (2,*) i,z(i)
 877 !        enddo
 878 !        call MPI_Barrier(FG_COMM,IERROR)
 879         time00=MPI_Wtime()
 880 !elwrite(iout,*) "do tej pory jest OK, MPI_Scatterv w ginv_mult"
 881         call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
 882           MPI_DOUBLE_PRECISION,&
 883           z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 884 !        write (2,*) "My chunk of z"
 885 !        do i=1,3*my_ng_count
 886 !          write (2,*) i,z(i)
 887 !        enddo
 888 !        write (2,*) "After SCATTERV"
 889 !        call flush(2)
 890 !        write (2,*) "MPI_Wtime",MPI_Wtime()
 891         time_scatter=time_scatter+MPI_Wtime()-time00
 892 #ifdef TIMING
 893         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
 894 #endif
 895 !        write (2,*) "time_scatter",time_scatter
 896 !        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
 897 !     &    my_ng_count
 898 !        call flush(2)
 899         time01=MPI_Wtime()
 900         do k=0,2
 901           do i=1,dimen
 902             ind=(i-1)*3+k+1
 903             temp(ind)=0.0d0
 904             do j=1,my_ng_count
 905 !              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
 906 !     &         Ginv(i,j),z((j-1)*3+k+1),
 907 !     &          Ginv(i,j)*z((j-1)*3+k+1)
 908 !              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
 909               temp(ind)=temp(ind)+Ginv(j,i)*z1((j-1)*3+k+1)
 910             enddo
 911           enddo
 912         enddo
 913         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 914 !        write (2,*) "Before REDUCE"
 915 !        call flush(2)
 916 !        write (2,*) "z before reduce"
 917 !        do i=1,dimen
 918 !          write (2,*) i,temp(i)
 919 !        enddo
 920         time00=MPI_Wtime()
 921         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
 922             MPI_SUM,king,FG_COMM,IERR)
 923         time_reduce=time_reduce+MPI_Wtime()-time00
 924 !        write (2,*) "After REDUCE"
 925 !        call flush(2)
 926       else
 927 #endif
 928 #ifdef TIMING
 929         time01=MPI_Wtime()
 930 #endif
 931         do k=0,2
 932           do i=1,dimen
 933             ind=(i-1)*3+k+1
 934             d_a_tmp(ind)=0.0d0
 935             do j=1,dimen
 936 !              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
 937 !              call flush(2)
 938 !     &         Ginv(i,j),z((j-1)*3+k+1),
 939 !     &          Ginv(i,j)*z((j-1)*3+k+1)
 940               d_a_tmp(ind)=d_a_tmp(ind) &
 941                                +Ginv(j,i)*z((j-1)*3+k+1)
 942 !              d_a_tmp(ind)=d_a_tmp(ind)
 943 !     &                         +Ginv(i,j)*z((j-1)*3+k+1)
 944             enddo
 945           enddo
 946         enddo
 947 #ifdef TIMING
 948         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 949 #endif
 950 #ifdef MPI
 951       endif
 952 #endif
 953       return
 954       end subroutine ginv_mult
 955 !-----------------------------------------------------------------------------
 956 #ifdef GINV_MULT
 957       subroutine ginv_mult_test(z,d_a_tmp)
 958
 959 !      include 'DIMENSIONS'
 960 !el      integer :: dimen
 961 !      include 'COMMON.MD'
 962       real(kind=8),dimension(dimen) :: z,d_a_tmp
 963       real(kind=8),dimension(dimen/3) :: ztmp,dtmp
 964       integer :: i,j,k,ind
 965 !      do i=1,dimen
 966 !        d_a_tmp(i)=0.0d0
 967 !        do j=1,dimen
 968 !          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
 969 !        enddo
 970 !      enddo
 971 !
 972 !      return
 973
 974 !ibm* unroll(3)
 975       do k=0,2
 976        do j=1,dimen/3
 977         ztmp(j)=z((j-1)*3+k+1)
 978        enddo
 979
 980        call alignx(16,ztmp(1))
 981        call alignx(16,dtmp(1))
 982        call alignx(16,Ginv(1,1))
 983
 984        do i=1,dimen/3
 985         dtmp(i)=0.0d0
 986         do j=1,dimen/3
 987            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
 988         enddo
 989        enddo
 990        do i=1,dimen/3
 991         ind=(i-1)*3+k+1
 992         d_a_tmp(ind)=dtmp(i)
 993        enddo
 994       enddo
 995       return
 996       end subroutine ginv_mult_test
 997 #endif
 998 !-----------------------------------------------------------------------------
 999       subroutine fricmat_mult(z,d_a_tmp)
1000
1001       use geometry_data, only: nres
1002       use control_data
1003       use MPI_data
1004 !      include 'DIMENSIONS'
1005 #ifdef MPI
1006       include 'mpif.h'
1007       integer :: IERROR,ierr
1008 #endif
1009 !      include 'COMMON.MD'
1010 !      include 'COMMON.IOUNITS'
1011 !      include 'COMMON.SETUP'
1012 !      include 'COMMON.TIME1'
1013 !#ifndef LANG0
1014 !      include 'COMMON.LANGEVIN'
1015 !#else
1016 !      include 'COMMON.LANGEVIN.lang0'
1017 !#endif
1018       real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
1019       real(kind=8),dimension(6*nres) :: temp    !(maxres6) maxres6=6*maxres
1020       real(kind=8) :: time00,time01
1021       integer :: i,j,k,ind,nres2
1022       nres2=2*nres
1023 !el      if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
1024
1025 #ifdef MPI
1026       if (nfgtasks.gt.1) then
1027         if (fg_rank.eq.0) then
1028 ! The matching BROADCAST for fg processors is called in ERGASTULUM
1029           time00=MPI_Wtime()
1030           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
1031           time_Bcast=time_Bcast+MPI_Wtime()-time00
1032 !          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
1033         endif
1034 !        call MPI_Barrier(FG_COMM,IERROR)
1035         time00=MPI_Wtime()
1036         call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
1037           MPI_DOUBLE_PRECISION,&
1038           z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1039 !        write (2,*) "My chunk of z"
1040 !        do i=1,3*my_ng_count
1041 !          write (2,*) i,z(i)
1042 !        enddo
1043         time_scatter=time_scatter+MPI_Wtime()-time00
1044 #ifdef TIMING
1045         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
1046 #endif
1047         time01=MPI_Wtime()
1048         do k=0,2
1049           do i=1,dimen
1050             ind=(i-1)*3+k+1
1051             temp(ind)=0.0d0
1052             do j=1,my_ng_count
1053               temp(ind)=temp(ind)-fricmat(j,i)*z1((j-1)*3+k+1)
1054             enddo
1055           enddo
1056         enddo
1057         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1058 !        write (2,*) "Before REDUCE"
1059 !        write (2,*) "d_a_tmp before reduce"
1060 !        do i=1,dimen3
1061 !          write (2,*) i,temp(i)
1062 !        enddo
1063 !        call flush(2)
1064         time00=MPI_Wtime()
1065         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
1066             MPI_SUM,king,FG_COMM,IERR)
1067         time_reduce=time_reduce+MPI_Wtime()-time00
1068 !        write (2,*) "After REDUCE"
1069 !        call flush(2)
1070       else
1071 #endif
1072 #ifdef TIMING
1073         time01=MPI_Wtime()
1074 #endif
1075         do k=0,2
1076          do i=1,dimen
1077           ind=(i-1)*3+k+1
1078           d_a_tmp(ind)=0.0d0
1079           do j=1,dimen
1080              d_a_tmp(ind)=d_a_tmp(ind) &
1081                                  -fricmat(j,i)*z((j-1)*3+k+1)
1082           enddo
1083          enddo
1084         enddo
1085 #ifdef TIMING
1086         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1087 #endif
1088 #ifdef MPI
1089       endif
1090 #endif
1091 !      write (iout,*) "Vector d_a"
1092 !      do i=1,dimen3
1093 !        write (2,*) i,d_a_tmp(i)
1094 !      enddo
1095       return
1096       end subroutine fricmat_mult
1097 !-----------------------------------------------------------------------------
1098 !-----------------------------------------------------------------------------
1099       end module REMD