source/unres/src_MD/thread.F

   1       subroutine thread_seq
   2 C Thread the sequence through a database of known structures
   3       implicit real*8 (a-h,o-z)
   4       include 'DIMENSIONS'
   5 #ifdef MPI
   6       include 'mpif.h'
   7 #endif
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.CHAIN'
  10       include 'COMMON.DBASE'
  11       include 'COMMON.INTERACT'
  12       include 'COMMON.VAR'
  13       include 'COMMON.THREAD'
  14       include 'COMMON.FFIELD'
  15       include 'COMMON.SBRIDGE'
  16       include 'COMMON.HEADER'
  17       include 'COMMON.IOUNITS'
  18       include 'COMMON.TIME1'
  19       include 'COMMON.CONTACTS'
  20       include 'COMMON.MCM'
  21       include 'COMMON.NAMES'
  22 #ifdef MPI
  23       include 'COMMON.INFO'
  24       integer ThreadId,ThreadType,Kwita
  25 #endif
  26       double precision varia(maxvar)
  27       double precision przes(3),obr(3,3)
  28       double precision time_for_thread
  29       logical found_pattern,non_conv
  30       character*32 head_pdb
  31       double precision energia(0:n_ene)
  32       n_ene_comp=nprint_ene
  33 C
  34 C Body
  35 C
  36 #ifdef MPI
  37       if (me.eq.king) then
  38         do i=1,nctasks
  39           nsave_part(i)=0
  40         enddo
  41       endif
  42       nacc_tot=0
  43 #endif
  44       Kwita=0
  45       close(igeom)
  46       close(ipdb)
  47       close(istat)
  48       do i=1,maxthread
  49         do j=1,14
  50           ener0(j,i)=0.0D0
  51           ener(j,i)=0.0D0
  52         enddo
  53       enddo
  54       nres0=nct-nnt+1
  55       ave_time_for_thread=0.0D0
  56       max_time_for_thread=0.0D0
  57 cd    print *,'nthread=',nthread,' nseq=',nseq,' nres0=',nres0
  58       nthread=nexcl+nthread
  59       do ithread=1,nthread
  60         found_pattern=.false.
  61         itrial=0
  62         do while (.not.found_pattern)
  63           itrial=itrial+1
  64           if (itrial.gt.1000) then
  65             write (iout,'(/a/)') 'Too many attempts to find pattern.'
  66             nthread=ithread-1
  67 #ifdef MPI
  68             call recv_stop_sig(Kwita)
  69             call send_stop_sig(-3)
  70 #endif
  71             goto 777
  72           endif
  73 C Find long enough chain in the database
  74           ii=iran_num(1,nseq)
  75           nres_t=nres_base(1,ii)
  76 C Select the starting position to thread.
  77           print *,'nseq',nseq,' ii=',ii,' nres_t=',
  78      &      nres_t,' nres0=',nres0
  79           if (nres_t.ge.nres0) then
  80             ist=iran_num(0,nres_t-nres0)
  81 #ifdef MPI
  82             if (Kwita.eq.0) call recv_stop_sig(Kwita)
  83             if (Kwita.lt.0) then
  84               write (iout,*) 'Stop signal received. Terminating.'
  85               write (*,*) 'Stop signal received. Terminating.'
  86               nthread=ithread-1
  87               write (*,*) 'ithread=',ithread,' nthread=',nthread
  88               goto 777
  89             endif
  90             call pattern_receive
  91 #endif
  92             do i=1,nexcl
  93               if (iexam(1,i).eq.ii .and. iexam(2,i).eq.ist) goto 10
  94             enddo
  95             found_pattern=.true.
  96           endif
  97 C If this point is reached, the pattern has not yet been examined.
  98    10     continue
  99 c         print *,'found_pattern:',found_pattern
 100         enddo
 101         nexcl=nexcl+1
 102         iexam(1,nexcl)=ii
 103         iexam(2,nexcl)=ist
 104 #ifdef MPI
 105         if (Kwita.eq.0) call recv_stop_sig(Kwita)
 106         if (Kwita.lt.0) then
 107           write (iout,*) 'Stop signal received. Terminating.'
 108           nthread=ithread-1
 109           write (*,*) 'ithread=',ithread,' nthread=',nthread
 110           goto 777
 111         endif
 112         call pattern_send
 113 #endif
 114         ipatt(1,ithread)=ii
 115         ipatt(2,ithread)=ist
 116 #ifdef MPI
 117         write (iout,'(/80(1h*)/a,i4,a,i5,2a,i3,a,i3,a,i3/)')
 118      &   'Processor:',me,' Attempt:',ithread,
 119      &   ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),
 120      &   ' start at res.',ist+1
 121         write (*,'(a,i4,a,i5,2a,i3,a,i3,a,i3)') 'Processor:',me,
 122      &   ' Attempt:',ithread,
 123      &   ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),
 124      &   ' start at res.',ist+1
 125 #else
 126         write (iout,'(/80(1h*)/a,i5,2a,i3,a,i3,a,i3/)')
 127      &   'Attempt:',ithread,
 128      &   ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),
 129      &   ' start at res.',ist+1
 130         write (*,'(a,i5,2a,i3,a,i3,a,i3)')
 131      &   'Attempt:',ithread,
 132      &   ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),
 133      &   ' start at res.',ist+1
 134 #endif
 135         ipattern=ii
 136 C Copy coordinates from the database.
 137         ist=ist-(nnt-1)
 138         do i=nnt,nct
 139           do j=1,3
 140             c(j,i)=cart_base(j,i+ist,ii)
 141 c           cref(j,i)=c(j,i)
 142           enddo
 143 cd        write (iout,'(a,i4,3f10.5)') restyp(itype(i)),i,(c(j,i),j=1,3)
 144         enddo
 145 cd      call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,
 146 cd             non_conv)
 147 cd      write (iout,'(a,f10.5)')
 148 cd   &  'Initial RMS deviation from reference structure:',rms
 149         if (itype(nres).eq.21) then
 150           do j=1,3
 151             dcj=c(j,nres-2)-c(j,nres-3)
 152             c(j,nres)=c(j,nres-1)+dcj
 153             c(j,2*nres)=c(j,nres)
 154           enddo
 155         endif
 156         if (itype(1).eq.21) then
 157           do j=1,3
 158             dcj=c(j,4)-c(j,3)
 159             c(j,1)=c(j,2)-dcj
 160             c(j,nres+1)=c(j,1)
 161           enddo
 162         endif
 163         call int_from_cart(.false.,.false.)
 164 cd      print *,'Exit INT_FROM_CART.'
 165 cd      print *,'nhpb=',nhpb
 166         do i=nss+1,nhpb
 167           ii=ihpb(i)
 168           jj=jhpb(i)
 169           dhpb(i)=dist(ii,jj)
 170 c         write (iout,'(2i5,2f10.5)') ihpb(i),jhpb(i),dhpb(i),forcon(i)
 171         enddo
 172 c       stop 'End generate'
 173 C Generate SC conformations.
 174         call sc_conf
 175 c       call intout
 176 #ifdef MPI
 177 cd      print *,'Processor:',me,': exit GEN_SIDE.'
 178 #else
 179 cd      print *,'Exit GEN_SIDE.'
 180 #endif
 181 C Calculate initial energy.
 182         call chainbuild
 183         call etotal(energia(0))
 184         etot=energia(0)
 185         do i=1,n_ene_comp
 186           ener0(i,ithread)=energia(i)
 187         enddo
 188         ener0(n_ene_comp+1,ithread)=energia(0)
 189         if (refstr) then
 190           call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
 191           ener0(n_ene_comp+3,ithread)=contact_fract(ncont,ncont_ref,
 192      &        icont,icont_ref)
 193           ener0(n_ene_comp+2,ithread)=rms
 194           ener0(n_ene_comp+4,ithread)=frac
 195           ener0(n_ene_comp+5,ithread)=frac_nn
 196         endif
 197         ener0(n_ene_comp+3,ithread)=0.0d0
 198 C Minimize energy.
 199 #ifdef MPI
 200        print*,'Processor:',me,' ithread=',ithread,' Start REGULARIZE.'
 201 #else
 202         print*,'ithread=',ithread,' Start REGULARIZE.'
 203 #endif
 204         curr_tim=tcpu()
 205         call regularize(nct-nnt+1,etot,rms,
 206      &                  cart_base(1,ist+nnt,ipattern),iretcode)
 207         curr_tim1=tcpu()
 208         time_for_thread=curr_tim1-curr_tim
 209         ave_time_for_thread=
 210      &  ((ithread-1)*ave_time_for_thread+time_for_thread)/ithread
 211         if (time_for_thread.gt.max_time_for_thread)
 212      &   max_time_for_thread=time_for_thread
 213 #ifdef MPI
 214         print *,'Processor',me,': Exit REGULARIZE.'
 215         if (WhatsUp.eq.2) then
 216           write (iout,*)
 217      &  'Sufficient number of confs. collected. Terminating.'
 218           nthread=ithread-1
 219           goto 777
 220         else if (WhatsUp.eq.-1) then
 221           nthread=ithread-1
 222           write (iout,*) 'Time up in REGULARIZE. Call SEND_STOP_SIG.'
 223           if (Kwita.eq.0) call recv_stop_sig(Kwita)
 224           call send_stop_sig(-2)
 225           goto 777
 226         else if (WhatsUp.eq.-2) then
 227           nthread=ithread-1
 228           write (iout,*) 'Timeup signal received. Terminating.'
 229           goto 777
 230         else if (WhatsUp.eq.-3) then
 231           nthread=ithread-1
 232           write (iout,*) 'Error stop signal received. Terminating.'
 233           goto 777
 234         endif
 235 #else
 236         print *,'Exit REGULARIZE.'
 237         if (iretcode.eq.11) then
 238           write (iout,'(/a/)')
 239      &'******* Allocated time exceeded in SUMSL. The program will stop.'
 240           nthread=ithread-1
 241           goto 777
 242         endif
 243 #endif
 244         head_pdb=titel(:24)//':'//str_nam(ipattern)
 245         if (outpdb) call pdbout(etot,head_pdb,ipdb)
 246         if (outmol2) call mol2out(etot,head_pdb)
 247 c       call intout
 248         call briefout(ithread,etot)
 249         link_end0=link_end
 250         link_end=min0(link_end,nss)
 251         write (iout,*) 'link_end=',link_end,' link_end0=',link_end0,
 252      &                 ' nss=',nss
 253         call etotal(energia(0))
 254 c       call enerprint(energia(0))
 255         link_end=link_end0
 256 cd      call chainbuild
 257 cd      call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,non_conv)
 258 cd      write (iout,'(a,f10.5)')
 259 cd   &  'RMS deviation from reference structure:',dsqrt(rms)
 260         do i=1,n_ene_comp
 261           ener(i,ithread)=energia(i)
 262         enddo
 263         ener(n_ene_comp+1,ithread)=energia(0)
 264         ener(n_ene_comp+3,ithread)=rms
 265         if (refstr) then
 266           call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
 267           ener(n_ene_comp+2,ithread)=rms
 268           ener(n_ene_comp+4,ithread)=frac
 269           ener(n_ene_comp+5,ithread)=frac_nn
 270         endif
 271         call write_stat_thread(ithread,ipattern,ist)
 272 c        write (istat,'(i4,2x,a8,i4,11(1pe14.5),2(0pf8.3),f8.5)')
 273 c     &  ithread,str_nam(ipattern),ist+1,(ener(k,ithread),k=1,11),
 274 c     &  (ener(k,ithread),k=12,14)
 275 #ifdef MPI
 276         if (me.eq.king) then
 277           nacc_tot=nacc_tot+1
 278           call pattern_receive
 279           call receive_MCM_info
 280           if (nacc_tot.ge.nthread) then
 281             write (iout,*)
 282      &     'Sufficient number of conformations collected nacc_tot=',
 283      &     nacc_tot,'. Stopping other processors and terminating.'
 284             write (*,*)
 285      &     'Sufficient number of conformations collected nacc_tot=',
 286      &     nacc_tot,'. Stopping other processors and terminating.'
 287            call recv_stop_sig(Kwita)
 288            if (Kwita.eq.0) call send_stop_sig(-1)
 289            nthread=ithread
 290            goto 777
 291           endif
 292         else
 293           call send_MCM_info(2)
 294         endif
 295 #endif
 296         if (timlim-curr_tim1-safety .lt. max_time_for_thread) then
 297           write (iout,'(/2a)')
 298      & '********** There would be not enough time for another thread. ',
 299      & 'The program will stop.'
 300           write (*,'(/2a)')
 301      & '********** There would be not enough time for another thread. ',
 302      & 'The program will stop.'
 303           write (iout,'(a,1pe14.4/)')
 304      &    'Elapsed time for last threading step: ',time_for_thread
 305           nthread=ithread
 306 #ifdef MPI
 307           call recv_stop_sig(Kwita)
 308           call send_stop_sig(-2)
 309 #endif
 310           goto 777
 311         else
 312           curr_tim=curr_tim1
 313           write (iout,'(a,1pe14.4)')
 314      &    'Elapsed time for this threading step: ',time_for_thread
 315         endif
 316 #ifdef MPI
 317         if (Kwita.eq.0) call recv_stop_sig(Kwita)
 318         if (Kwita.lt.0) then
 319           write (iout,*) 'Stop signal received. Terminating.'
 320           write (*,*) 'Stop signal received. Terminating.'
 321           nthread=ithread
 322           write (*,*) 'nthread=',nthread,' ithread=',ithread
 323           goto 777
 324         endif
 325 #endif
 326       enddo
 327 #ifdef MPI
 328       call send_stop_sig(-1)
 329 #endif
 330   777 continue
 331 #ifdef MPI
 332 C Any messages left for me?
 333       call pattern_receive
 334       if (Kwita.eq.0) call recv_stop_sig(Kwita)
 335 #endif
 336       call write_thread_summary
 337 #ifdef MPI
 338       if (king.eq.king) then
 339         Kwita=1
 340         do while (Kwita.ne.0 .or. nacc_tot.ne.0)
 341           Kwita=0
 342           nacc_tot=0
 343           call recv_stop_sig(Kwita)
 344           call receive_MCM_info
 345         enddo
 346         do iproc=1,nprocs-1
 347           call receive_thread_results(iproc)
 348         enddo
 349         call write_thread_summary
 350       else
 351         call send_thread_results
 352       endif
 353 #endif
 354       return
 355       end
 356 c--------------------------------------------------------------------------
 357       subroutine write_thread_summary
 358 C Thread the sequence through a database of known structures
 359       implicit real*8 (a-h,o-z)
 360       include 'DIMENSIONS'
 361 #ifdef MPI
 362       include 'mpif.h'
 363 #endif
 364       include 'COMMON.CONTROL'
 365       include 'COMMON.CHAIN'
 366       include 'COMMON.DBASE'
 367       include 'COMMON.INTERACT'
 368       include 'COMMON.VAR'
 369       include 'COMMON.THREAD'
 370       include 'COMMON.FFIELD'
 371       include 'COMMON.SBRIDGE'
 372       include 'COMMON.HEADER'
 373       include 'COMMON.NAMES'
 374       include 'COMMON.IOUNITS'
 375       include 'COMMON.TIME1'
 376 #ifdef MPI
 377       include 'COMMON.INFO'
 378 #endif
 379       dimension ip(maxthread)
 380       double precision energia(0:n_ene)
 381       write (iout,'(30x,a/)')
 382      & '  *********** Summary threading statistics ************'
 383       write (iout,'(a)') 'Initial energies:'
 384       write (iout,'(a4,2x,a12,14a14,3a8)')
 385      & 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',
 386      & 'RMSnat','NatCONT','NNCONT','RMS'
 387 C Energy sort patterns
 388       do i=1,nthread
 389         ip(i)=i
 390       enddo
 391       do i=1,nthread-1
 392         enet=ener(n_ene-1,ip(i))
 393         jj=i
 394         do j=i+1,nthread
 395           if (ener(n_ene-1,ip(j)).lt.enet) then
 396             jj=j
 397             enet=ener(n_ene-1,ip(j))
 398           endif
 399         enddo
 400         if (jj.ne.i) then
 401           ipj=ip(jj)
 402           ip(jj)=ip(i)
 403           ip(i)=ipj
 404         endif
 405       enddo
 406       do ik=1,nthread
 407         i=ip(ik)
 408         ii=ipatt(1,i)
 409         ist=nres_base(2,ii)+ipatt(2,i)
 410         do kk=1,n_ene_comp
 411           energia(i)=ener0(kk,i)
 412         enddo
 413         etot=ener0(n_ene_comp+1,i)
 414         rmsnat=ener0(n_ene_comp+2,i)
 415         rms=ener0(n_ene_comp+3,i)
 416         frac=ener0(n_ene_comp+4,i)
 417         frac_nn=ener0(n_ene_comp+5,i)
 418
 419         if (refstr) then
 420         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)')
 421      &  i,str_nam(ii),ist+1,
 422      &  (energia(print_order(kk)),kk=1,nprint_ene),
 423      &  etot,rmsnat,frac,frac_nn,rms
 424         else
 425         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)')
 426      &  i,str_nam(ii),ist+1,
 427      &  (energia(print_order(kk)),kk=1,nprint_ene),etot
 428         endif
 429       enddo
 430       write (iout,'(//a)') 'Final energies:'
 431       write (iout,'(a4,2x,a12,17a14,3a8)')
 432      & 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',
 433      & 'RMSnat','NatCONT','NNCONT','RMS'
 434       do ik=1,nthread
 435         i=ip(ik)
 436         ii=ipatt(1,i)
 437         ist=nres_base(2,ii)+ipatt(2,i)
 438         do kk=1,n_ene_comp
 439           energia(kk)=ener(kk,ik)
 440         enddo
 441         etot=ener(n_ene_comp+1,i)
 442         rmsnat=ener(n_ene_comp+2,i)
 443         rms=ener(n_ene_comp+3,i)
 444         frac=ener(n_ene_comp+4,i)
 445         frac_nn=ener(n_ene_comp+5,i)
 446         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)')
 447      &  i,str_nam(ii),ist+1,
 448      &  (energia(print_order(kk)),kk=1,nprint_ene),
 449      &  etot,rmsnat,frac,frac_nn,rms
 450       enddo
 451       write (iout,'(/a/)') 'IEXAM array:'
 452       write (iout,'(i5)') nexcl
 453       do i=1,nexcl
 454         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
 455       enddo
 456       write (iout,'(/a,1pe14.4/a,1pe14.4/)')
 457      & 'Max. time for threading step ',max_time_for_thread,
 458      & 'Average time for threading step: ',ave_time_for_thread
 459       return
 460       end
 461 c----------------------------------------------------------------------------
 462       subroutine sc_conf
 463 C Sample (hopefully) optimal SC orientations given backcone conformation.
 464       implicit real*8 (a-h,o-z)
 465       include 'DIMENSIONS'
 466       include 'COMMON.CHAIN'
 467       include 'COMMON.DBASE'
 468       include 'COMMON.INTERACT'
 469       include 'COMMON.VAR'
 470       include 'COMMON.THREAD'
 471       include 'COMMON.FFIELD'
 472       include 'COMMON.SBRIDGE'
 473       include 'COMMON.HEADER'
 474       include 'COMMON.GEO'
 475       include 'COMMON.IOUNITS'
 476       double precision varia(maxvar)
 477       common /srutu/ icall
 478       double precision energia(0:n_ene)
 479       logical glycine,fail
 480       maxsample=10
 481       link_end0=link_end
 482       link_end=min0(link_end,nss)
 483       do i=nnt,nct
 484         if (itype(i).ne.10) then
 485 cd        print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1)
 486           call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail)
 487         endif
 488       enddo
 489       call chainbuild
 490       call etotal(energia(0))
 491       do isample=1,maxsample
 492 C Choose a non-glycine side chain.
 493         glycine=.true.
 494         do while(glycine)
 495           ind_sc=iran_num(nnt,nct)
 496           glycine=(itype(ind_sc).eq.10)
 497         enddo
 498         alph0=alph(ind_sc)
 499         omeg0=omeg(ind_sc)
 500         call gen_side(itype(ind_sc),theta(ind_sc+1),alph(ind_sc),
 501      &       omeg(ind_sc),fail)
 502         call chainbuild
 503         call etotal(energia(0))
 504 cd      write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))')
 505 cd   &   'Step:',isample,' SC',ind_sc,' alpha',alph(ind_sc)*rad2deg,
 506 cd   &   ' omega',omeg(ind_sc)*rad2deg,' old energy',e0,' new energy',e1
 507         e1=0.0d0
 508         if (e0.le.e1) then
 509           alph(ind_sc)=alph0
 510           omeg(ind_sc)=omeg0
 511         else
 512           e0=e1
 513         endif
 514       enddo
 515       link_end=link_end0
 516       return
 517       end
 518 c---------------------------------------------------------------------------
 519       subroutine write_stat_thread(ithread,ipattern,ist)
 520       implicit real*8 (a-h,o-z)
 521       include "DIMENSIONS"
 522       include "COMMON.CONTROL"
 523       include "COMMON.IOUNITS"
 524       include "COMMON.THREAD"
 525       include "COMMON.FFIELD"
 526       include "COMMON.DBASE"
 527       include "COMMON.NAMES"
 528       double precision energia(0:n_ene)
 529
 530 #if defined(AIX) || defined(PGI)
 531       open(istat,file=statname,position='append')
 532 #else
 533       open(istat,file=statname,access='append')
 534 #endif
 535       do i=1,n_ene_comp
 536         energia(i)=ener(i,ithread)
 537       enddo
 538       etot=ener(n_ene_comp+1,ithread)
 539       rmsnat=ener(n_ene_comp+2,ithread)
 540       rms=ener(n_ene_comp+3,ithread)
 541       frac=ener(n_ene_comp+4,ithread)
 542       frac_nn=ener(n_ene_comp+5,ithread)
 543       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)')
 544      &  ithread,str_nam(ipattern),ist+1,
 545      &  (energia(print_order(i)),i=1,nprint_ene),
 546      &  etot,rmsnat,frac,frac_nn,rms
 547       close (istat)
 548       return
 549       end