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