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