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 = .false.
 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       print *,"just entered"
 323       relfeh=1.0d-14
 324       nres2=2*nres
 325       print *,"FIVEDIAG"
 326       write (iout,*) "before FIVEDIAG"
 327 #ifndef FIVEDIAG
 328       write (iout,*) "ALLOCATE"
 329       print *,"ALLOCATE"
 330       if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
 331       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
 332 #endif
 333 !
 334 ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
 335 ! inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
 336 !
 337 ! Determine the number of degrees of freedom (dimen) and the number of
 338 ! sites (dimen1)
 339       dimen=(nct-nnt+1)+nside
 340       dimen1=(nct-nnt)+(nct-nnt+1)
 341       dimen3=dimen*3
 342       write (iout,*) "nnt",nnt," nct",nct," nside",nside
 343 #ifdef FIVEDIAG
 344        ip4=ip/4
 345        DM(1)=mp/4+ip4
 346        if (iabs(itype(nnt)).eq.10) then
 347          DM(1)=DM(1)+msc(10)
 348          ind=2
 349          nind=1
 350        else
 351          DM(1)=DM(1)+isc(iabs(itype(nnt)))
 352          DM(2)=msc(iabs(itype(nnt)))+isc(iabs(itype(nnt)))
 353          ind=3
 354          nind=2
 355        endif
 356        do i=nnt+1,nct-1
 357 !         if (iabs(itype(i,1)).eq.ntyp1) cycle
 358          DM(ind)=2*ip4+mp/2
 359          if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) then
 360            if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10)
 361            ind=ind+1
 362          else
 363            DM(ind)=DM(ind)+isc(iabs(itype(i,1)))
 364            DM(ind+1)=msc(iabs(itype(i,1)))+isc(iabs(itype(i,1)))
 365            ind=ind+2
 366          endif
 367        enddo
 368        if (nct.gt.nnt) then
 369        DM(ind)=ip4+mp/4
 370        if (iabs(itype(nct)).eq.10) then
 371          DM(ind)=DM(ind)+msc(10)
 372          nind=ind
 373        else
 374          DM(ind)=DM(ind)+isc(iabs(itype(nct)))
 375          DM(ind+1)=msc(iabs(itype(nct)))+isc(iabs(itype(nct)))
 376          nind=ind+1
 377         endif
 378         endif
 379
 380
 381         ind=1
 382         do i=nnt,nct
 383         mnum=molnum(i)
 384         if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 &
 385           .and.mnum.eq.5) then
 386           DU1(ind)=-isc(iabs(itype(i,1)))
 387           DU1(ind+1)=0.0d0
 388           ind=ind+2
 389         else
 390           DU1(ind)=mp/4-ip4
 391           ind=ind+1
 392         endif
 393         enddo
 394
 395        ind=1
 396        do i=nnt,nct-1
 397        mnum=molnum(i)
 398 !        if (iabs(itype(i,1)).eq.ntyp1) cycle
 399         write (iout,*) "i",i," itype",itype(i,1),ntyp1
 400         if (iabs(itype(i,1)).ne.10 .and. &
 401           iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
 402          DU2(ind)=mp/4-ip4
 403          DU2(ind+1)=0.0d0
 404          ind=ind+2
 405         else
 406          DU2(ind)=0.0d0
 407          DU2(ind+1)=0.0d0
 408          ind=ind+1
 409         endif
 410        enddo
 411        DMorig=DM
 412        DU1orig=DU1
 413        DU2orig=DU2
 414        write (iout,*) "nind",nind," dimen",dimen
 415        if (nind.ne.dimen) then
 416          write (iout,*) "ERROR, dimensions differ"
 417 #ifdef MPI
 418          call MPI_Finalize(ierr)
 419 #endif
 420          stop
 421        endif
 422      write (iout,*) "The upper part of the five-diagonal inertia matrix"
 423        do i=1,dimen
 424          if (i.lt.dimen-1) then
 425            write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
 426          else if (i.eq.dimen-1) then
 427            write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
 428          else
 429            write (iout,'(i3,3f10.5)') i,DM(i)
 430          endif
 431        enddo
 432
 433         call FDISYP (dimen, DM, DU1, DU2, MARK)
 434
 435         if (mark.eq.-1) then
 436      write (iout,*) "ERROR: the inertia matrix is not positive definite"
 437 #ifdef MPI
 438          call MPI_Finalize(ierr)
 439 #endif
 440          stop
 441         else if (mark.eq.0) then
 442      write (iout,*) "ERROR: the inertia matrix is singular"
 443 #ifdef MPI
 444          call MPI_Finalize(ierr)
 445 #endif
 446         else if (mark.eq.1) then
 447           write (iout,*) "The transformed inertia matrix"
 448           do i=1,dimen
 449             if (i.lt.dimen-1) then
 450               write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
 451             else if (i.eq.dimen-1) then
 452               write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
 453             else
 454               write (iout,'(i3,3f10.5)') i,DM(i)
 455             endif
 456           enddo
 457         endif
 458 ! Diagonalization of the pentadiagonal matrix
 459         allocate(DDU1(2*nres))
 460         allocate(DDU2(2*nres))
 461         allocate(DDM(2*nres))
 462         do i=1,dimen-1
 463           DDU1(i+1)=DU1orig(i)
 464         enddo
 465         do i=1,dimen-2
 466           DDU2(i+2)=DU2orig(i)
 467         enddo
 468         DDM=DMorig
 469         call quindibisect2(DDM,DDU1,DDU2,dimen,1,dimen,eps1,relfeh,eps2,iz,geigen)
 470         if (lprn) then
 471           write (iout,*) &
 472            "Eigenvalues of the five-diagonal inertia matrix"
 473           do i=1,dimen
 474             write (iout,'(i5,f10.5)') i,geigen(i)
 475           enddo
 476         endif
 477         deallocate(DDU1)
 478         deallocate(DDU2)
 479         deallocate(DDM)
 480 #else
 481 ! Old Gmatrix
 482 #ifdef MPI
 483       if (nfgtasks.gt.1) then
 484       time00=MPI_Wtime()
 485       call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
 486       time_Bcast=time_Bcast+MPI_Wtime()-time00
 487       call int_bounds(dimen,igmult_start,igmult_end)
 488       igmult_start=igmult_start-1
 489       call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,&
 490           ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 491       my_ng_count=igmult_end-igmult_start
 492       call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,&
 493           MPI_INTEGER,FG_COMM,IERROR)
 494       write (iout,*) 'Processor:',fg_rank,' CG group',kolor,&
 495        ' absolute rank',myrank,' igmult_start',igmult_start,&
 496        ' igmult_end',igmult_end,' count',my_ng_count
 497       write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 498       write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 499       call flush(iout)
 500       else
 501 #endif
 502       igmult_start=1
 503       igmult_end=dimen
 504       my_ng_count=dimen
 505 #ifdef MPI
 506       endif
 507 #endif
 508 !      write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
 509 !  Zeroing out A and fricmat
 510       do i=1,dimen
 511         do j=1,dimen
 512           A(i,j)=0.0D0
 513         enddo
 514       enddo
 515 !  Diagonal elements of the dC part of A and the respective friction coefficients
 516       ind=1
 517       ind1=0
 518       print *,"TUTUTUT?!",nnt,nct-1
 519       do i=nnt,nct-1
 520         mnum=molnum(i)
 521         ind=ind+1
 522         ind1=ind1+1
 523         coeff=0.25d0*IP(mnum)
 524         if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 525         print *,i,coeff,mp(mnum)
 526         massvec(ind1)=mp(mnum)
 527         Gmat(ind,ind)=coeff
 528         print *,"i",mp(mnum)
 529         A(ind1,ind)=0.5d0
 530       enddo
 531
 532 !  Off-diagonal elements of the dC part of A
 533       k=3
 534       do i=1,nct-nnt
 535         do j=1,i
 536           A(i,j)=1.0d0
 537         enddo
 538       enddo
 539 !  Diagonal elements of the dX part of A and the respective friction coefficients
 540       m=nct-nnt
 541       m1=nct-nnt+1
 542       ind=0
 543       ind1=0
 544       do i=1,2
 545       msc(ntyp1_molec(i),i)=1.0d0
 546       enddo
 547       do i=nnt,nct
 548         ind=ind+1
 549         ii = ind+m
 550         mnum=molnum(i)
 551         iti=itype(i,mnum)
 552         if (mnum.eq.5) then
 553         mscab=0.0
 554         else
 555         mscab=msc(iabs(iti),mnum)
 556         endif
 557         massvec(ii)=mscab
 558
 559         if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
 560           ind1=ind1+1
 561           ii1= ind1+m1
 562           A(ii,ii1)=1.0d0
 563           Gmat(ii1,ii1)=ISC(iabs(iti),mnum)
 564         endif
 565       enddo
 566 !  Off-diagonal elements of the dX part of A
 567       ind=0
 568       k=nct-nnt
 569       do i=nnt,nct
 570         iti=itype(i,1)
 571         ind=ind+1
 572         do j=nnt,i
 573           ii = ind
 574           jj = j-nnt+1
 575           A(k+ii,jj)=1.0d0
 576         enddo
 577       enddo
 578       if (lprn) then
 579         write (iout,*)
 580         write (iout,*) "Vector massvec"
 581         do i=1,dimen1
 582           write (iout,*) i,massvec(i)
 583         enddo
 584         write (iout,'(//a)') "A"
 585         call matout(dimen,dimen1,nres2,nres2,A)
 586       endif
 587
 588 ! Calculate the G matrix (store in Gmat)
 589       do k=1,dimen
 590        do i=1,dimen
 591          dtdi=0.0d0
 592          do j=1,dimen1
 593            dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
 594          enddo
 595          Gmat(k,i)=Gmat(k,i)+dtdi
 596        enddo
 597       enddo
 598
 599       if (lprn) then
 600         write (iout,'(//a)') "Gmat"
 601         call matout(dimen,dimen,nres2,nres2,Gmat)
 602       endif
 603       do i=1,dimen
 604         do j=1,dimen
 605           Ginv(i,j)=0.0d0
 606           Gcopy(i,j)=Gmat(i,j)
 607         enddo
 608         Ginv(i,i)=1.0d0
 609       enddo
 610 ! Invert the G matrix
 611       call MATINVERT(dimen,nres2,Gcopy,Ginv,osob)
 612       if (lprn) then
 613         write (iout,'(//a)') "Ginv"
 614         call matout(dimen,dimen,nres2,nres2,Ginv)
 615       endif
 616
 617       Bmat=0.0d0
 618       Bmat(1,1)=1.0d0
 619       j=2
 620       do i=nnt,nct
 621         mnum=molnum(i)
 622         if ((itype(i,1).eq.10).or.(itype(i,mnum).eq.ntyp1_molec(mnum))&
 623           .or.(mnum.eq.5)) then
 624           Bmat(i-nnt+2,j-1)=-1
 625           Bmat(i-nnt+2,j)=1
 626           j=j+1
 627         else
 628           if (i.lt.nct) then
 629             Bmat(i-nnt+2,j-1)=-1
 630             Bmat(i-nnt+2,j+1)=1
 631             Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
 632             Bmat(i-nnt+1+(nct-nnt)+1,j)=1
 633             j=j+2
 634           else
 635             Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
 636             Bmat(i-nnt+1+(nct-nnt)+1,j)=1
 637           endif
 638         endif
 639       enddo
 640       write (iout,*) "j",j," dimen",dimen
 641       if (lprn) then
 642         write (iout,'(//a)') "Bmat"
 643         call matout(dimen,dimen,nres2,nres2,Bmat)
 644       endif
 645       Gcopy=0.0d0
 646       do i=1,dimen
 647         do j=1,dimen
 648             do k=1,dimen
 649                do l=1,dimen
 650                  Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j)
 651                enddo
 652             enddo
 653         enddo
 654       enddo
 655       if (lprn) then
 656         write (iout,'(//a)') "Gmat-transformed"
 657         call matout(dimen,dimen,nres2,nres2,Gcopy)
 658       endif
 659 #ifdef MPI
 660       if (nfgtasks.gt.1) then
 661         myginv_ng_count=nres2*my_ng_count
 662         call MPI_Allgather(nres2*igmult_start,1,MPI_INTEGER,&
 663           nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 664         call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,&
 665           nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 666         write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 667         write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 668         call flush(iout)
 669 !        call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 670 !     &    nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
 671 !     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 672 !        call MPI_Barrier(FG_COMM,IERR)
 673         time00=MPI_Wtime()
 674         call MPI_Scatterv(ginv(1,1),nginv_counts(0),&
 675           nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),&
 676           myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 677 #ifdef TIMING
 678         time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 679 #endif
 680         do i=1,dimen
 681           do j=1,2*my_ng_count
 682             ginv(j,i)=gcopy(i,j)
 683           enddo
 684         enddo
 685 !        write (iout,*) "Master's chunk of ginv"
 686 !        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
 687       endif
 688 #endif
 689       if (osob) then
 690         write (iout,*) "The G matrix is singular."
 691         stop
 692       endif
 693 ! Compute G**(-1/2) and G**(1/2)
 694       ind=0
 695       do i=1,dimen
 696         do j=1,i
 697           ind=ind+1
 698           Ghalf(ind)=Gmat(i,j)
 699         enddo
 700       enddo
 701       call gldiag(nres2,dimen,dimen,Ghalf,work,Geigen,Gvec,&
 702         ierr,iwork)
 703       if (lprn) then
 704         write (iout,'(//a)') &
 705          "Eigenvectors and eigenvalues of the G matrix"
 706         call eigout(dimen,dimen,nres2,nres2,Gvec,Geigen)
 707       endif
 708       do i=1,dimen
 709         sqreig(i)=dsqrt(Geigen(i))
 710       enddo
 711       do i=1,dimen
 712         do j=1,dimen
 713           Gsqrp(i,j)=0.0d0
 714           Gsqrm(i,j)=0.0d0
 715           Gcopy(i,j)=0.0d0
 716           do k=1,dimen
 717             Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
 718             Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
 719             Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
 720           enddo
 721         enddo
 722       enddo
 723       if (lprn) then
 724         write (iout,*) "Comparison of original and restored G"
 725         do i=1,dimen
 726           do j=1,dimen
 727             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),&
 728               Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
 729           enddo
 730         enddo
 731       endif
 732       deallocate(Gcopy)
 733 #endif
 734       return
 735       end subroutine setup_MD_matrices
 736 !-----------------------------------------------------------------------------
 737       subroutine EIGOUT(NC,NR,LM2,LM3,A,B)
 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),B(LM2)
 744       KA=1
 745       KC=6
 746     1 KB=MIN0(KC,NC)
 747       WRITE(IOUT,600) (I,I=KA,KB)
 748       WRITE(IOUT,601) (B(I),I=KA,KB)
 749       WRITE(IOUT,602)
 750     2 N=0
 751       DO 3  I=1,NR
 752       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 753       N=N+1
 754       IF(N.LT.10) GO TO 3
 755       WRITE(IOUT,602)
 756       N=0
 757     3 CONTINUE
 758     4 IF (KB.EQ.NC) RETURN
 759       KA=KC+1
 760       KC=KC+6
 761       GO TO 1
 762   600 FORMAT (// 9H ROOT NO.,I4,9I11)
 763   601 FORMAT (/5X,10(1PE11.4))
 764   602 FORMAT (2H  )
 765   603 FORMAT (I5,10F11.5)
 766   604 FORMAT (1H1)
 767       end subroutine EIGOUT
 768 !-----------------------------------------------------------------------------
 769       subroutine MATOUT(NC,NR,LM2,LM3,A)
 770
 771 !      implicit real*8 (a-h,o-z)
 772 !      include 'DIMENSIONS'
 773 !      include 'COMMON.IOUNITS'
 774       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 775       real(kind=8) :: A(LM2,LM3)
 776       KA=1
 777       KC=6
 778     1 KB=MIN0(KC,NC)
 779       WRITE(IOUT,600) (I,I=KA,KB)
 780       WRITE(IOUT,602)
 781     2 N=0
 782       DO 3  I=1,NR
 783       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 784       N=N+1
 785       IF(N.LT.10) GO TO 3
 786       WRITE(IOUT,602)
 787       N=0
 788     3 CONTINUE
 789     4 IF (KB.EQ.NC) RETURN
 790       KA=KC+1
 791       KC=KC+6
 792       GO TO 1
 793   600 FORMAT (//5x,9I11)
 794   602 FORMAT (2H  )
 795   603 FORMAT (I5,10F11.3)
 796   604 FORMAT (1H1)
 797       end subroutine MATOUT
 798 !-----------------------------------------------------------------------------
 799       subroutine MATOUT1(NC,NR,LM2,LM3,A)
 800
 801 !      implicit real*8 (a-h,o-z)
 802 !      include 'DIMENSIONS'
 803 !      include 'COMMON.IOUNITS'
 804       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 805       real(kind=8) :: A(LM2,LM3)
 806       KA=1
 807       KC=21
 808     1 KB=MIN0(KC,NC)
 809       WRITE(IOUT,600) (I,I=KA,KB)
 810       WRITE(IOUT,602)
 811     2 N=0
 812       DO 3  I=1,NR
 813       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 814       N=N+1
 815       IF(N.LT.3) GO TO 3
 816       WRITE(IOUT,602)
 817       N=0
 818     3 CONTINUE
 819     4 IF (KB.EQ.NC) RETURN
 820       KA=KC+1
 821       KC=KC+21
 822       GO TO 1
 823   600 FORMAT (//5x,7(3I5,2x))
 824   602 FORMAT (2H  )
 825   603 FORMAT (I5,7(3F5.1,2x))
 826   604 FORMAT (1H1)
 827       end subroutine MATOUT1
 828 !-----------------------------------------------------------------------------
 829       subroutine MATOUT2(NC,NR,LM2,LM3,A)
 830
 831 !      implicit real*8 (a-h,o-z)
 832 !      include 'DIMENSIONS'
 833 !      include 'COMMON.IOUNITS'
 834       integer :: I,J,KA,KC,KB,N
 835       integer :: LM2,LM3,NC,NR
 836       real(kind=8) :: A(LM2,LM3)
 837       KA=1
 838       KC=12
 839     1 KB=MIN0(KC,NC)
 840       WRITE(IOUT,600) (I,I=KA,KB)
 841       WRITE(IOUT,602)
 842     2 N=0
 843       DO 3  I=1,NR
 844       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 845       N=N+1
 846       IF(N.LT.3) GO TO 3
 847       WRITE(IOUT,602)
 848       N=0
 849     3 CONTINUE
 850     4 IF (KB.EQ.NC) RETURN
 851       KA=KC+1
 852       KC=KC+12
 853       GO TO 1
 854   600 FORMAT (//5x,4(3I9,2x))
 855   602 FORMAT (2H  )
 856   603 FORMAT (I5,4(3F9.3,2x))
 857   604 FORMAT (1H1)
 858       end subroutine MATOUT2
 859 !-----------------------------------------------------------------------------
 860       subroutine ginv_mult(z,d_a_tmp)
 861
 862       use geometry_data, only: nres
 863       use control_data
 864       use MPI_data
 865 !      implicit real*8 (a-h,o-z)
 866 !      include 'DIMENSIONS'
 867 #ifdef MPI
 868       include 'mpif.h'
 869       integer :: ierr,ierror
 870 #endif
 871 !      include 'COMMON.SETUP'
 872 !      include 'COMMON.TIME1'
 873 !      include 'COMMON.MD'
 874       real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
 875       real(kind=8),dimension(6*nres) :: temp    !(maxres6) maxres6=6*maxres
 876       real(kind=8) :: time00,time01
 877       integer :: i,j,k,ind
 878 #ifdef MPI
 879       if (nfgtasks.gt.1) then
 880         if (fg_rank.eq.0) then
 881 ! The matching BROADCAST for fg processors is called in ERGASTULUM
 882           time00=MPI_Wtime()
 883           call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
 884           time_Bcast=time_Bcast+MPI_Wtime()-time00
 885 !          print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
 886         endif
 887 !        write (2,*) "time00",time00
 888 !        write (2,*) "Before Scatterv"
 889 !        call flush(2)
 890 !        write (2,*) "Whole z (for FG master)"
 891 !        do i=1,dimen
 892 !          write (2,*) i,z(i)
 893 !        enddo
 894 !        call MPI_Barrier(FG_COMM,IERROR)
 895         time00=MPI_Wtime()
 896 !elwrite(iout,*) "do tej pory jest OK, MPI_Scatterv w ginv_mult"
 897         call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
 898           MPI_DOUBLE_PRECISION,&
 899           z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 900 !        write (2,*) "My chunk of z"
 901 !        do i=1,3*my_ng_count
 902 !          write (2,*) i,z(i)
 903 !        enddo
 904 !        write (2,*) "After SCATTERV"
 905 !        call flush(2)
 906 !        write (2,*) "MPI_Wtime",MPI_Wtime()
 907         time_scatter=time_scatter+MPI_Wtime()-time00
 908 #ifdef TIMING
 909         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
 910 #endif
 911 !        write (2,*) "time_scatter",time_scatter
 912 !        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
 913 !     &    my_ng_count
 914 !        call flush(2)
 915         time01=MPI_Wtime()
 916         do k=0,2
 917           do i=1,dimen
 918             ind=(i-1)*3+k+1
 919             temp(ind)=0.0d0
 920             do j=1,my_ng_count
 921 !              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
 922 !     &         Ginv(i,j),z((j-1)*3+k+1),
 923 !     &          Ginv(i,j)*z((j-1)*3+k+1)
 924 !              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
 925               temp(ind)=temp(ind)+Ginv(j,i)*z1((j-1)*3+k+1)
 926             enddo
 927           enddo
 928         enddo
 929         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 930 !        write (2,*) "Before REDUCE"
 931 !        call flush(2)
 932 !        write (2,*) "z before reduce"
 933 !        do i=1,dimen
 934 !          write (2,*) i,temp(i)
 935 !        enddo
 936         time00=MPI_Wtime()
 937         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
 938             MPI_SUM,king,FG_COMM,IERR)
 939         time_reduce=time_reduce+MPI_Wtime()-time00
 940 !        write (2,*) "After REDUCE"
 941 !        call flush(2)
 942       else
 943 #endif
 944 #ifdef TIMING
 945         time01=MPI_Wtime()
 946 #endif
 947         do k=0,2
 948           do i=1,dimen
 949             ind=(i-1)*3+k+1
 950             d_a_tmp(ind)=0.0d0
 951             do j=1,dimen
 952 !              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
 953 !              call flush(2)
 954 !     &         Ginv(i,j),z((j-1)*3+k+1),
 955 !     &          Ginv(i,j)*z((j-1)*3+k+1)
 956               d_a_tmp(ind)=d_a_tmp(ind) &
 957                                +Ginv(j,i)*z((j-1)*3+k+1)
 958 !              d_a_tmp(ind)=d_a_tmp(ind)
 959 !     &                         +Ginv(i,j)*z((j-1)*3+k+1)
 960             enddo
 961           enddo
 962         enddo
 963 #ifdef TIMING
 964         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 965 #endif
 966 #ifdef MPI
 967       endif
 968 #endif
 969       return
 970       end subroutine ginv_mult
 971 !-----------------------------------------------------------------------------
 972 #ifdef GINV_MULT
 973       subroutine ginv_mult_test(z,d_a_tmp)
 974
 975 !      include 'DIMENSIONS'
 976 !el      integer :: dimen
 977 !      include 'COMMON.MD'
 978       real(kind=8),dimension(dimen) :: z,d_a_tmp
 979       real(kind=8),dimension(dimen/3) :: ztmp,dtmp
 980       integer :: i,j,k,ind
 981 !      do i=1,dimen
 982 !        d_a_tmp(i)=0.0d0
 983 !        do j=1,dimen
 984 !          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
 985 !        enddo
 986 !      enddo
 987 !
 988 !      return
 989
 990 !ibm* unroll(3)
 991       do k=0,2
 992        do j=1,dimen/3
 993         ztmp(j)=z((j-1)*3+k+1)
 994        enddo
 995
 996        call alignx(16,ztmp(1))
 997        call alignx(16,dtmp(1))
 998        call alignx(16,Ginv(1,1))
 999
1000        do i=1,dimen/3
1001         dtmp(i)=0.0d0
1002         do j=1,dimen/3
1003            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
1004         enddo
1005        enddo
1006        do i=1,dimen/3
1007         ind=(i-1)*3+k+1
1008         d_a_tmp(ind)=dtmp(i)
1009        enddo
1010       enddo
1011       return
1012       end subroutine ginv_mult_test
1013 #endif
1014 !-----------------------------------------------------------------------------
1015       subroutine fricmat_mult(z,d_a_tmp)
1016
1017       use geometry_data, only: nres
1018       use control_data
1019       use MPI_data
1020 !      include 'DIMENSIONS'
1021 #ifdef MPI
1022       include 'mpif.h'
1023       integer :: IERROR,ierr
1024 #endif
1025 !      include 'COMMON.MD'
1026 !      include 'COMMON.IOUNITS'
1027 !      include 'COMMON.SETUP'
1028 !      include 'COMMON.TIME1'
1029 !#ifndef LANG0
1030 !      include 'COMMON.LANGEVIN'
1031 !#else
1032 !      include 'COMMON.LANGEVIN.lang0'
1033 !#endif
1034       real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
1035       real(kind=8),dimension(6*nres) :: temp    !(maxres6) maxres6=6*maxres
1036       real(kind=8) :: time00,time01
1037       integer :: i,j,k,ind,nres2
1038       nres2=2*nres
1039 !el      if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
1040
1041 #ifdef MPI
1042       if (nfgtasks.gt.1) then
1043         if (fg_rank.eq.0) then
1044 ! The matching BROADCAST for fg processors is called in ERGASTULUM
1045           time00=MPI_Wtime()
1046           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
1047           time_Bcast=time_Bcast+MPI_Wtime()-time00
1048 !          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
1049         endif
1050 !        call MPI_Barrier(FG_COMM,IERROR)
1051         time00=MPI_Wtime()
1052         call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
1053           MPI_DOUBLE_PRECISION,&
1054           z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1055 !        write (2,*) "My chunk of z"
1056 !        do i=1,3*my_ng_count
1057 !          write (2,*) i,z(i)
1058 !        enddo
1059         time_scatter=time_scatter+MPI_Wtime()-time00
1060 #ifdef TIMING
1061         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
1062 #endif
1063         time01=MPI_Wtime()
1064         do k=0,2
1065           do i=1,dimen
1066             ind=(i-1)*3+k+1
1067             temp(ind)=0.0d0
1068             do j=1,my_ng_count
1069               temp(ind)=temp(ind)-fricmat(j,i)*z1((j-1)*3+k+1)
1070             enddo
1071           enddo
1072         enddo
1073         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1074 !        write (2,*) "Before REDUCE"
1075 !        write (2,*) "d_a_tmp before reduce"
1076 !        do i=1,dimen3
1077 !          write (2,*) i,temp(i)
1078 !        enddo
1079 !        call flush(2)
1080         time00=MPI_Wtime()
1081         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
1082             MPI_SUM,king,FG_COMM,IERR)
1083         time_reduce=time_reduce+MPI_Wtime()-time00
1084 !        write (2,*) "After REDUCE"
1085 !        call flush(2)
1086       else
1087 #endif
1088 #ifdef TIMING
1089         time01=MPI_Wtime()
1090 #endif
1091         do k=0,2
1092          do i=1,dimen
1093           ind=(i-1)*3+k+1
1094           d_a_tmp(ind)=0.0d0
1095           do j=1,dimen
1096              d_a_tmp(ind)=d_a_tmp(ind) &
1097                                  -fricmat(j,i)*z((j-1)*3+k+1)
1098           enddo
1099          enddo
1100         enddo
1101 #ifdef TIMING
1102         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1103 #endif
1104 #ifdef MPI
1105       endif
1106 #endif
1107 !      write (iout,*) "Vector d_a"
1108 !      do i=1,dimen3
1109 !        write (2,*) i,d_a_tmp(i)
1110 !      enddo
1111       return
1112       end subroutine fricmat_mult
1113 !-----------------------------------------------------------------------------
1114 !-----------------------------------------------------------------------------
1115       end module REMD