source/unres/src_CSA/unres_csa.F

   1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
   2 C                                                                              C
   3 C                                U N R E S                                     C
   4 C                                                                              C
   5 C Program to carry out conformational search of proteins in an united-residue  C
   6 C approximation.                                                               C
   7 C                                                                              C
   8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
   9       implicit real*8 (a-h,o-z)
  10       include 'DIMENSIONS'
  11
  12
  13 #ifdef MPI
  14       include 'mpif.h'
  15       include 'COMMON.SETUP'
  16 #endif
  17       include 'COMMON.TIME1'
  18       include 'COMMON.INTERACT'
  19       include 'COMMON.NAMES'
  20       include 'COMMON.GEO'
  21       include 'COMMON.HEADER'
  22       include 'COMMON.CONTROL'
  23       include 'COMMON.CONTACTS'
  24       include 'COMMON.CHAIN'
  25       include 'COMMON.VAR'
  26       include 'COMMON.IOUNITS'
  27       include 'COMMON.FFIELD'
  28 c      include 'COMMON.REMD'
  29 c      include 'COMMON.MD'
  30       include 'COMMON.SBRIDGE'
  31       double precision hrtime,mintime,sectime
  32       character*64 text_mode_calc(-2:14) /'test',
  33      & 'SC rotamer distribution',
  34      & 'Energy evaluation or minimization',
  35      & 'Regularization of PDB structure',
  36      & 'Threading of a sequence on PDB structures',
  37      & 'Monte Carlo (with minimization) ',
  38      & 'Energy minimization of multiple conformations',
  39      & 'Checking energy gradient',
  40      & 'Entropic sampling Monte Carlo (with minimization)',
  41      & 'Energy map',
  42      & 'CSA calculations',
  43      & 'Not used 9',
  44      & 'Not used 10',
  45      & 'Soft regularization of PDB structure',
  46      & 'Mesoscopic molecular dynamics (MD) ',
  47      & 'Not used 13',
  48      & 'Replica exchange molecular dynamics (REMD)'/
  49       external ilen
  50
  51 c      call memmon_print_usage()
  52
  53       call init_task
  54       if (me.eq.king)
  55      & write(iout,*)'### LAST MODIFIED  11/03/09 1:19PM by czarek'
  56       if (me.eq.king) call cinfo
  57 C Read force field parameters and job setup data
  58       call readrtns
  59       call flush(iout)
  60 C
  61       if (me.eq.king .or. .not. out1file) then
  62        write (iout,'(2a/)')
  63      & text_mode_calc(modecalc)(:ilen(text_mode_calc(modecalc))),
  64      & ' calculation.'
  65        if (minim) write (iout,'(a)')
  66      &  'Conformations will be energy-minimized.'
  67        write (iout,'(80(1h*)/)')
  68       endif
  69       call flush(iout)
  70 C
  71 c      if (modecalc.eq.-2) then
  72 c        call test
  73 c        stop
  74 c      else if (modecalc.eq.-1) then
  75 c        write(iout,*) "call check_sc_map next"
  76 c        call check_bond
  77 c        stop
  78 c      endif
  79 #ifdef MPI
  80       if (fg_rank.gt.0) then
  81 C Fine-grain slaves just do energy and gradient components.
  82         call ergastulum ! slave workhouse in Latin
  83       else
  84 #endif
  85       if (modecalc.eq.0) then
  86         call exec_eeval_or_minim
  87 c      else if (modecalc.eq.1) then
  88 c        call exec_regularize
  89 c      else if (modecalc.eq.2) then
  90 c        call exec_thread
  91 c      else if (modecalc.eq.3 .or. modecalc .eq.6) then
  92 c        call exec_MC
  93       else if (modecalc.eq.4) then
  94         call exec_mult_eeval_or_minim
  95       else if (modecalc.eq.5) then
  96          call exec_checkgrad
  97 c      else if (ModeCalc.eq.7) then
  98 c        call exec_map
  99       else if (ModeCalc.eq.8) then
 100         call exec_CSA
 101 c      else if (modecalc.eq.11) then
 102 c        call exec_softreg
 103 c      else if (modecalc.eq.12) then
 104 c        call exec_MD
 105 c      else if (modecalc.eq.14) then
 106 c        call exec_MREMD
 107       else
 108         write (iout,'(a)') 'This calculation type is not supported',
 109      &   ModeCalc
 110       endif
 111 #ifdef MPI
 112       endif
 113 C Finish task.
 114       if (fg_rank.eq.0) call finish_task
 115 c      call memmon_print_usage()
 116 #ifdef TIMING
 117        call print_detailed_timing
 118 #endif
 119       call MPI_Finalize(ierr)
 120       stop 'Bye Bye...'
 121 #else
 122       call dajczas(tcpu(),hrtime,mintime,sectime)
 123       stop '********** Program terminated normally.'
 124 #endif
 125       end
 126 c---------------------------------------------------------------------------
 127       subroutine exec_eeval_or_minim
 128       implicit real*8 (a-h,o-z)
 129       include 'DIMENSIONS'
 130 #ifdef MPI
 131       include 'mpif.h'
 132 #endif
 133       include 'COMMON.SETUP'
 134       include 'COMMON.TIME1'
 135       include 'COMMON.INTERACT'
 136       include 'COMMON.NAMES'
 137       include 'COMMON.GEO'
 138       include 'COMMON.HEADER'
 139       include 'COMMON.CONTROL'
 140       include 'COMMON.CONTACTS'
 141       include 'COMMON.CHAIN'
 142       include 'COMMON.VAR'
 143       include 'COMMON.IOUNITS'
 144       include 'COMMON.FFIELD'
 145 c      include 'COMMON.REMD'
 146 c      include 'COMMON.MD'
 147       include 'COMMON.SBRIDGE'
 148       common /srutu/ icall
 149       double precision energy(0:n_ene)
 150       double precision energy_long(0:n_ene),energy_short(0:n_ene)
 151       double precision varia(maxvar)
 152       if (indpdb.eq.0) call chainbuild
 153 c      time00=MPI_Wtime()
 154       call chainbuild_cart
 155 c      if (split_ene) then
 156 c       print *,"Processor",myrank," after chainbuild"
 157 c       icall=1
 158 c       call etotal_long(energy_long(0))
 159 c       write (iout,*) "Printing long range energy"
 160 c       call enerprint(energy_long(0))
 161 c       call etotal_short(energy_short(0))
 162 c       write (iout,*) "Printing short range energy"
 163 c       call enerprint(energy_short(0))
 164 c       do i=0,n_ene
 165 c         energy(i)=energy_long(i)+energy_short(i)
 166 c         write (iout,*) i,energy_long(i),energy_short(i),energy(i)
 167 c       enddo
 168 c       write (iout,*) "Printing long+short range energy"
 169 c       call enerprint(energy(0))
 170 c      endif
 171       call etotal(energy(0))
 172 c      time_ene=MPI_Wtime()-time00
 173       write (iout,*) "Time for energy evaluation",time_ene
 174       print *,"after etotal"
 175       etota = energy(0)
 176       etot =etota
 177       call enerprint(energy(0))
 178 c      call hairpin(.true.,nharp,iharp)
 179 c      call secondary2(.true.)
 180       if (minim) then
 181
 182         if (dccart) then
 183           print *, 'Calling MINIM_DC'
 184 c          time1=MPI_WTIME()
 185           call minim_dc(etot,iretcode,nfun)
 186         else
 187           if (indpdb.ne.0) then
 188             call bond_regular
 189             call chainbuild
 190           endif
 191           call geom_to_var(nvar,varia)
 192           print *,'Calling MINIMIZE.'
 193 c          time1=MPI_WTIME()
 194           call minimize(etot,varia,iretcode,nfun)
 195         endif
 196         print *,'SUMSL return code is',iretcode,' eval ',nfun
 197 c        evals=nfun/(MPI_WTIME()-time1)
 198         print *,'# eval/s',evals
 199         print *,'refstr=',refstr
 200 c        call hairpin(.true.,nharp,iharp)
 201 c        call secondary2(.true.)
 202         call etotal(energy(0))
 203         etot = energy(0)
 204         call enerprint(energy(0))
 205
 206         call intout
 207         call briefout(0,etot)
 208         if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
 209           write (iout,'(a,i3)') 'SUMSL return code:',iretcode
 210           write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
 211           write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
 212       else
 213         print *,'refstr=',refstr
 214         if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
 215         call briefout(0,etot)
 216       endif
 217       if (outpdb) call pdbout(etot,titel(:32),ipdb)
 218       if (outmol2) call mol2out(etot,titel(:32))
 219       return
 220       end
 221
 222       subroutine exec_checkgrad
 223       implicit real*8 (a-h,o-z)
 224       include 'DIMENSIONS'
 225 #ifdef MPI
 226       include 'mpif.h'
 227 #endif
 228       include 'COMMON.SETUP'
 229       include 'COMMON.TIME1'
 230       include 'COMMON.INTERACT'
 231       include 'COMMON.NAMES'
 232       include 'COMMON.GEO'
 233       include 'COMMON.HEADER'
 234       include 'COMMON.CONTROL'
 235       include 'COMMON.CONTACTS'
 236       include 'COMMON.CHAIN'
 237       include 'COMMON.VAR'
 238       include 'COMMON.IOUNITS'
 239       include 'COMMON.FFIELD'
 240 c      include 'COMMON.REMD'
 241       include 'COMMON.MD_'
 242       include 'COMMON.SBRIDGE'
 243       common /srutu/ icall
 244       double precision energy(0:max_ene)
 245 c      do i=2,nres
 246 c        vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0)
 247 c        if (itype(i).ne.10)
 248 c     &      vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0)
 249 c      enddo
 250       if (indpdb.eq.0) call chainbuild
 251 c      do i=0,nres
 252 c        do j=1,3
 253 c          dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0)
 254 c        enddo
 255 c      enddo
 256 c      do i=1,nres-1
 257 c        if (itype(i).ne.10) then
 258 c          do j=1,3
 259 c            dc(j,i+nres)=dc(j,i+nres)+ran_number(-0.2d0,0.2d0)
 260 c          enddo
 261 c        endif
 262 c      enddo
 263 c      do j=1,3
 264 c        dc(j,0)=ran_number(-0.2d0,0.2d0)
 265 c      enddo
 266       usampl=.true.
 267       totT=1.d0
 268       eq_time=0.0d0
 269 c      call read_fragments
 270 c      read(inp,*) t_bath
 271 c      call rescale_weights(t_bath)
 272       call chainbuild_cart
 273       call cartprint
 274       call intout
 275       icall=1
 276       call etotal(energy(0))
 277       etot = energy(0)
 278       call enerprint(energy(0))
 279       write (iout,*) "Uconst",Uconst," Uconst_back",uconst_back
 280       print *,'icheckgrad=',icheckgrad
 281       goto (10,20,30) icheckgrad
 282   10  call check_ecartint
 283       return
 284   20  call check_cartgrad
 285       return
 286   30  call check_eint
 287       return
 288       end
 289 c---------------------------------------------------------------------------
 290       subroutine exec_CSA
 291 #ifdef MPI
 292       include "mpif.h"
 293 #endif
 294       include 'DIMENSIONS'
 295       include 'COMMON.IOUNITS'
 296 C Conformational Space Annealling programmed by Jooyoung Lee.
 297 C This method works only with parallel machines!
 298 #ifdef MPI
 299       call together
 300 #else
 301       write (iout,*) "CSA works on parallel machines only"
 302 #endif
 303       return
 304       end
 305 c---------------------------------------------------------------------------
 306 #ifdef MPI
 307       subroutine exec_mult_eeval_or_minim
 308       implicit real*8 (a-h,o-z)
 309       include 'DIMENSIONS'
 310       include 'mpif.h'
 311       integer muster(mpi_status_size)
 312       include 'COMMON.SETUP'
 313       include 'COMMON.TIME1'
 314       include 'COMMON.INTERACT'
 315       include 'COMMON.NAMES'
 316       include 'COMMON.GEO'
 317       include 'COMMON.HEADER'
 318       include 'COMMON.CONTROL'
 319       include 'COMMON.CONTACTS'
 320       include 'COMMON.CHAIN'
 321       include 'COMMON.VAR'
 322       include 'COMMON.IOUNITS'
 323       include 'COMMON.FFIELD'
 324       include 'COMMON.SBRIDGE'
 325       double precision varia(maxvar)
 326       integer ind(6)
 327       double precision energy(0:n_ene)
 328       logical eof
 329       eof=.false.
 330
 331       if(me.ne.king) then
 332         call minim_mcmf
 333         return
 334       endif
 335
 336       close (intin)
 337       open(intin,file=intinname,status='old')
 338       write (istat,'(a5,100a12)')"#    ",
 339      &  (wname(print_order(i)),i=1,nprint_ene)
 340       if (refstr) then
 341         write (istat,'(a5,100a12)')"#    ",
 342      &   (ename(print_order(i)),i=1,nprint_ene),
 343      &   "ETOT total","RMSD","nat.contact","nnt.contact",
 344      &   "cont.order","TMscore"
 345       else
 346         write (istat,'(a5,100a12)')"#    ",
 347      &    (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
 348       endif
 349
 350       if (.not.minim) then
 351         do while (.not. eof)
 352           if (read_cart) then
 353             read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene
 354             call read_x(intin,*11)
 355 c Broadcast the order to compute internal coordinates to the slaves.
 356             if (nfgtasks.gt.1)
 357      &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
 358             call int_from_cart1(.false.)
 359           else
 360             read (intin,'(i5)',end=1100,err=1100) iconf
 361             call read_angles(intin,*11)
 362             call geom_to_var(nvar,varia)
 363             call chainbuild
 364           endif
 365           write (iout,'(a,i7)') 'Conformation #',iconf
 366           call etotal(energy(0))
 367           call briefout(iconf,energy(0))
 368           call enerprint(energy(0))
 369           etot=energy(0)
 370           if (refstr) then
 371             call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
 372             call calc_tmscore(tm,.true.)
 373             write (istat,'(i5,100(f12.3))') iconf,
 374      &      (energy(print_order(i)),i=1,nprint_ene),etot,
 375      &       rms,frac,frac_nn,co,tm
 376           else
 377             write (istat,'(i5,100(f12.3))') iconf,
 378      &     (energy(print_order(i)),i=1,nprint_ene),etot
 379           endif
 380         enddo
 381 1100    continue
 382         goto 1101
 383       endif
 384
 385       mm=0
 386       imm=0
 387       nft=0
 388       ene0=0.0d0
 389       n=0
 390       iconf=0
 391       do while (.not. eof)
 392         mm=mm+1
 393         if (mm.lt.nodes) then
 394           if (read_cart) then
 395             read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
 396             call read_x(intin,*11)
 397 c Broadcast the order to compute internal coordinates to the slaves.
 398             if (nfgtasks.gt.1)
 399      &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
 400             call int_from_cart1(.false.)
 401           else
 402             read (intin,'(i5)',end=11,err=11) iconf
 403             call read_angles(intin,*11)
 404             call geom_to_var(nvar,varia)
 405             call chainbuild
 406           endif
 407
 408           n=n+1
 409           write (iout,*) 'Conformation #',iconf,' read'
 410          imm=imm+1
 411          ind(1)=1
 412          ind(2)=n
 413          ind(3)=0
 414          ind(4)=0
 415          ind(5)=0
 416          ind(6)=0
 417          ene0=0.0d0
 418          call mpi_send(ind,6,mpi_integer,mm,idint,CG_COMM,
 419      *                  ierr)
 420          call mpi_send(varia,nvar,mpi_double_precision,mm,
 421      *                  idreal,CG_COMM,ierr)
 422          call mpi_send(ene0,1,mpi_double_precision,mm,
 423      *                  idreal,CG_COMM,ierr)
 424 c         print *,'task ',n,' sent to worker ',mm,nvar
 425         else
 426          call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
 427      *                 CG_COMM,muster,ierr)
 428          man=muster(mpi_source)
 429 c         print *,'receiving result from worker ',man,' (',iii1,iii,')'
 430          call mpi_recv(varia,nvar,mpi_double_precision,
 431      *               man,idreal,CG_COMM,muster,ierr)
 432          call mpi_recv(ene,1,
 433      *               mpi_double_precision,man,idreal,
 434      *               CG_COMM,muster,ierr)
 435          call mpi_recv(ene0,1,
 436      *               mpi_double_precision,man,idreal,
 437      *               CG_COMM,muster,ierr)
 438 c         print *,'result received from worker ',man,' sending now'
 439
 440           call var_to_geom(nvar,varia)
 441           call chainbuild
 442           call etotal(energy(0))
 443           iconf=ind(2)
 444           write (iout,*)
 445           write (iout,*)
 446           write (iout,*) 'Conformation #',iconf,ind(5)
 447
 448           etot=energy(0)
 449           call enerprint(energy(0))
 450           call briefout(iconf,etot)
 451           if (refstr) then
 452             call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
 453             call calc_tmscore(tm,.true.)
 454             write (istat,'(i5,100(f12.3))') iconf,
 455      &     (energy(print_order(i)),i=1,nprint_ene),etot,
 456      &     rms,frac,frac_nn,co,tm
 457           else
 458             write (istat,'(i5,100(f12.3))') iconf,
 459      &     (energy(print_order(i)),i=1,nprint_ene),etot
 460           endif
 461
 462           imm=imm-1
 463           if (read_cart) then
 464             read (intin,'(e15.10,e15.5)',end=1101,err=1101) time,ene
 465             call read_x(intin,*11)
 466 c Broadcast the order to compute internal coordinates to the slaves.
 467             if (nfgtasks.gt.1)
 468      &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
 469             call int_from_cart1(.false.)
 470           else
 471             read (intin,'(i5)',end=1101,err=1101) iconf
 472             call read_angles(intin,*11)
 473             call geom_to_var(nvar,varia)
 474             call chainbuild
 475           endif
 476           n=n+1
 477           write (iout,*) 'Conformation #',iconf,' read'
 478           imm=imm+1
 479           ind(1)=1
 480           ind(2)=n
 481           ind(3)=0
 482           ind(4)=0
 483           ind(5)=0
 484           ind(6)=0
 485           call mpi_send(ind,6,mpi_integer,man,idint,CG_COMM,
 486      *                  ierr)
 487           call mpi_send(varia,nvar,mpi_double_precision,man,
 488      *                  idreal,CG_COMM,ierr)
 489           call mpi_send(ene0,1,mpi_double_precision,man,
 490      *                  idreal,CG_COMM,ierr)
 491           nf_mcmf=nf_mcmf+ind(4)
 492           nmin=nmin+1
 493         endif
 494       enddo
 495 11    continue
 496       do j=1,imm
 497         call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
 498      *               CG_COMM,muster,ierr)
 499         man=muster(mpi_source)
 500         call mpi_recv(varia,nvar,mpi_double_precision,
 501      *               man,idreal,CG_COMM,muster,ierr)
 502         call mpi_recv(ene,1,
 503      *               mpi_double_precision,man,idreal,
 504      *               CG_COMM,muster,ierr)
 505         call mpi_recv(ene0,1,
 506      *               mpi_double_precision,man,idreal,
 507      *               CG_COMM,muster,ierr)
 508
 509         call var_to_geom(nvar,varia)
 510         call chainbuild
 511         call etotal(energy(0))
 512         iconf=ind(2)
 513         write (iout,*)
 514         write (iout,*)
 515         write (iout,*) 'Conformation #',iconf,ind(5)
 516
 517         etot=energy(0)
 518         call enerprint(energy(0))
 519         call briefout(iconf,etot)
 520         if (refstr) then
 521           call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
 522           call calc_tmscore(tm,.true.)
 523           write (istat,'(i5,100(f12.3))') iconf,
 524      &   (energy(print_order(i)),i=1,nprint_ene),etot,
 525      &   rms,frac,frac_nn,co,tm
 526         else
 527           write (istat,'(i5,100(f12.3))') iconf,
 528      &    (energy(print_order(i)),i=1,nprint_ene),etot
 529         endif
 530         nmin=nmin+1
 531       enddo
 532 1101  continue
 533       do i=1, nodes-1
 534          ind(1)=0
 535          ind(2)=0
 536          ind(3)=0
 537          ind(4)=0
 538          ind(5)=0
 539          ind(6)=0
 540          call mpi_send(ind,6,mpi_integer,i,idint,CG_COMM,
 541      *                  ierr)
 542       enddo
 543       return
 544       end
 545 #else
 546       subroutine exec_mult_eeval_or_minim
 547       include 'DIMENSIONS'
 548       include 'COMMON.IOUNITS'
 549       write (iout,*) "Unsupported option in serial version"
 550       return
 551       end
 552 #endif
 553 c---------------------------------------------------------------------------
 554