--- /dev/null
+ subroutine thread_seq
+C Thread the sequence through a database of known structures
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DBASE'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.THREAD'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.HEADER'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.TIME1'
+ include 'COMMON.CONTACTS'
+ include 'COMMON.MCM'
+ include 'COMMON.NAMES'
+#ifdef MPI
+ include 'COMMON.INFO'
+ integer ThreadId,ThreadType,Kwita
+#endif
+ double precision varia(maxvar)
+ double precision przes(3),obr(3,3)
+ double precision time_for_thread
+ logical found_pattern,non_conv
+ character*32 head_pdb
+ double precision energia(0:n_ene)
+ n_ene_comp=nprint_ene
+C
+C Body
+C
+#ifdef MPI
+ if (me.eq.king) then
+ do i=1,nctasks
+ nsave_part(i)=0
+ enddo
+ endif
+ nacc_tot=0
+#endif
+ Kwita=0
+ close(igeom)
+ close(ipdb)
+ close(istat)
+ do i=1,maxthread
+ do j=1,14
+ ener0(j,i)=0.0D0
+ ener(j,i)=0.0D0
+ enddo
+ enddo
+ nres0=nct-nnt+1
+ ave_time_for_thread=0.0D0
+ max_time_for_thread=0.0D0
+cd print *,'nthread=',nthread,' nseq=',nseq,' nres0=',nres0
+ nthread=nexcl+nthread
+ do ithread=1,nthread
+ found_pattern=.false.
+ itrial=0
+ do while (.not.found_pattern)
+ itrial=itrial+1
+ if (itrial.gt.1000) then
+ write (iout,'(/a/)') 'Too many attempts to find pattern.'
+ nthread=ithread-1
+#ifdef MPI
+ call recv_stop_sig(Kwita)
+ call send_stop_sig(-3)
+#endif
+ goto 777
+ endif
+C Find long enough chain in the database
+ ii=iran_num(1,nseq)
+ nres_t=nres_base(1,ii)
+C Select the starting position to thread.
+ print *,'nseq',nseq,' ii=',ii,' nres_t=',
+ & nres_t,' nres0=',nres0
+ if (nres_t.ge.nres0) then
+ ist=iran_num(0,nres_t-nres0)
+#ifdef MPI
+ if (Kwita.eq.0) call recv_stop_sig(Kwita)
+ if (Kwita.lt.0) then
+ write (iout,*) 'Stop signal received. Terminating.'
+ write (*,*) 'Stop signal received. Terminating.'
+ nthread=ithread-1
+ write (*,*) 'ithread=',ithread,' nthread=',nthread
+ goto 777
+ endif
+ call pattern_receive
+#endif
+ do i=1,nexcl
+ if (iexam(1,i).eq.ii .and. iexam(2,i).eq.ist) goto 10
+ enddo
+ found_pattern=.true.
+ endif
+C If this point is reached, the pattern has not yet been examined.
+ 10 continue
+c print *,'found_pattern:',found_pattern
+ enddo
+ nexcl=nexcl+1
+ iexam(1,nexcl)=ii
+ iexam(2,nexcl)=ist
+#ifdef MPI
+ if (Kwita.eq.0) call recv_stop_sig(Kwita)
+ if (Kwita.lt.0) then
+ write (iout,*) 'Stop signal received. Terminating.'
+ nthread=ithread-1
+ write (*,*) 'ithread=',ithread,' nthread=',nthread
+ goto 777
+ endif
+ call pattern_send
+#endif
+ ipatt(1,ithread)=ii
+ ipatt(2,ithread)=ist
+#ifdef MPI
+ write (iout,'(/80(1h*)/a,i4,a,i5,2a,i3,a,i3,a,i3/)')
+ & 'Processor:',me,' Attempt:',ithread,
+ & ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),
+ & ' start at res.',ist+1
+ write (*,'(a,i4,a,i5,2a,i3,a,i3,a,i3)') 'Processor:',me,
+ & ' Attempt:',ithread,
+ & ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),
+ & ' start at res.',ist+1
+#else
+ write (iout,'(/80(1h*)/a,i5,2a,i3,a,i3,a,i3/)')
+ & 'Attempt:',ithread,
+ & ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),
+ & ' start at res.',ist+1
+ write (*,'(a,i5,2a,i3,a,i3,a,i3)')
+ & 'Attempt:',ithread,
+ & ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),
+ & ' start at res.',ist+1
+#endif
+ ipattern=ii
+C Copy coordinates from the database.
+ ist=ist-(nnt-1)
+ do i=nnt,nct
+ do j=1,3
+ c(j,i)=cart_base(j,i+ist,ii)
+c cref(j,i)=c(j,i)
+ enddo
+cd write (iout,'(a,i4,3f10.5)') restyp(itype(i)),i,(c(j,i),j=1,3)
+ enddo
+cd call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,
+cd non_conv)
+cd write (iout,'(a,f10.5)')
+cd & 'Initial RMS deviation from reference structure:',rms
+ if (itype(nres).eq.21) then
+ do j=1,3
+ dcj=c(j,nres-2)-c(j,nres-3)
+ c(j,nres)=c(j,nres-1)+dcj
+ c(j,2*nres)=c(j,nres)
+ enddo
+ endif
+ if (itype(1).eq.21) then
+ do j=1,3
+ dcj=c(j,4)-c(j,3)
+ c(j,1)=c(j,2)-dcj
+ c(j,nres+1)=c(j,1)
+ enddo
+ endif
+ call int_from_cart(.false.,.false.)
+cd print *,'Exit INT_FROM_CART.'
+cd print *,'nhpb=',nhpb
+ do i=nss+1,nhpb
+ ii=ihpb(i)
+ jj=jhpb(i)
+ dhpb(i)=dist(ii,jj)
+c write (iout,'(2i5,2f10.5)') ihpb(i),jhpb(i),dhpb(i),forcon(i)
+ enddo
+c stop 'End generate'
+C Generate SC conformations.
+ call sc_conf
+c call intout
+#ifdef MPI
+cd print *,'Processor:',me,': exit GEN_SIDE.'
+#else
+cd print *,'Exit GEN_SIDE.'
+#endif
+C Calculate initial energy.
+ call chainbuild
+ call etotal(energia(0))
+ etot=energia(0)
+ do i=1,n_ene_comp
+ ener0(i,ithread)=energia(i)
+ enddo
+ ener0(n_ene_comp+1,ithread)=energia(0)
+ if (refstr) then
+ call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+ ener0(n_ene_comp+3,ithread)=contact_fract(ncont,ncont_ref,
+ & icont,icont_ref)
+ ener0(n_ene_comp+2,ithread)=rms
+ ener0(n_ene_comp+4,ithread)=frac
+ ener0(n_ene_comp+5,ithread)=frac_nn
+ endif
+ ener0(n_ene_comp+3,ithread)=0.0d0
+C Minimize energy.
+#ifdef MPI
+ print*,'Processor:',me,' ithread=',ithread,' Start REGULARIZE.'
+#else
+ print*,'ithread=',ithread,' Start REGULARIZE.'
+#endif
+ curr_tim=tcpu()
+ call regularize(nct-nnt+1,etot,rms,
+ & cart_base(1,ist+nnt,ipattern),iretcode)
+ curr_tim1=tcpu()
+ time_for_thread=curr_tim1-curr_tim
+ ave_time_for_thread=
+ & ((ithread-1)*ave_time_for_thread+time_for_thread)/ithread
+ if (time_for_thread.gt.max_time_for_thread)
+ & max_time_for_thread=time_for_thread
+#ifdef MPI
+ print *,'Processor',me,': Exit REGULARIZE.'
+ if (WhatsUp.eq.2) then
+ write (iout,*)
+ & 'Sufficient number of confs. collected. Terminating.'
+ nthread=ithread-1
+ goto 777
+ else if (WhatsUp.eq.-1) then
+ nthread=ithread-1
+ write (iout,*) 'Time up in REGULARIZE. Call SEND_STOP_SIG.'
+ if (Kwita.eq.0) call recv_stop_sig(Kwita)
+ call send_stop_sig(-2)
+ goto 777
+ else if (WhatsUp.eq.-2) then
+ nthread=ithread-1
+ write (iout,*) 'Timeup signal received. Terminating.'
+ goto 777
+ else if (WhatsUp.eq.-3) then
+ nthread=ithread-1
+ write (iout,*) 'Error stop signal received. Terminating.'
+ goto 777
+ endif
+#else
+ print *,'Exit REGULARIZE.'
+ if (iretcode.eq.11) then
+ write (iout,'(/a/)')
+ &'******* Allocated time exceeded in SUMSL. The program will stop.'
+ nthread=ithread-1
+ goto 777
+ endif
+#endif
+ head_pdb=titel(:24)//':'//str_nam(ipattern)
+ if (outpdb) call pdbout(etot,head_pdb,ipdb)
+ if (outmol2) call mol2out(etot,head_pdb)
+c call intout
+ call briefout(ithread,etot)
+ link_end0=link_end
+ link_end=min0(link_end,nss)
+ write (iout,*) 'link_end=',link_end,' link_end0=',link_end0,
+ & ' nss=',nss
+ call etotal(energia(0))
+c call enerprint(energia(0))
+ link_end=link_end0
+cd call chainbuild
+cd call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,non_conv)
+cd write (iout,'(a,f10.5)')
+cd & 'RMS deviation from reference structure:',dsqrt(rms)
+ do i=1,n_ene_comp
+ ener(i,ithread)=energia(i)
+ enddo
+ ener(n_ene_comp+1,ithread)=energia(0)
+ ener(n_ene_comp+3,ithread)=rms
+ if (refstr) then
+ call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+ ener(n_ene_comp+2,ithread)=rms
+ ener(n_ene_comp+4,ithread)=frac
+ ener(n_ene_comp+5,ithread)=frac_nn
+ endif
+ call write_stat_thread(ithread,ipattern,ist)
+c write (istat,'(i4,2x,a8,i4,11(1pe14.5),2(0pf8.3),f8.5)')
+c & ithread,str_nam(ipattern),ist+1,(ener(k,ithread),k=1,11),
+c & (ener(k,ithread),k=12,14)
+#ifdef MPI
+ if (me.eq.king) then
+ nacc_tot=nacc_tot+1
+ call pattern_receive
+ call receive_MCM_info
+ if (nacc_tot.ge.nthread) then
+ write (iout,*)
+ & 'Sufficient number of conformations collected nacc_tot=',
+ & nacc_tot,'. Stopping other processors and terminating.'
+ write (*,*)
+ & 'Sufficient number of conformations collected nacc_tot=',
+ & nacc_tot,'. Stopping other processors and terminating.'
+ call recv_stop_sig(Kwita)
+ if (Kwita.eq.0) call send_stop_sig(-1)
+ nthread=ithread
+ goto 777
+ endif
+ else
+ call send_MCM_info(2)
+ endif
+#endif
+ if (timlim-curr_tim1-safety .lt. max_time_for_thread) then
+ write (iout,'(/2a)')
+ & '********** There would be not enough time for another thread. ',
+ & 'The program will stop.'
+ write (*,'(/2a)')
+ & '********** There would be not enough time for another thread. ',
+ & 'The program will stop.'
+ write (iout,'(a,1pe14.4/)')
+ & 'Elapsed time for last threading step: ',time_for_thread
+ nthread=ithread
+#ifdef MPI
+ call recv_stop_sig(Kwita)
+ call send_stop_sig(-2)
+#endif
+ goto 777
+ else
+ curr_tim=curr_tim1
+ write (iout,'(a,1pe14.4)')
+ & 'Elapsed time for this threading step: ',time_for_thread
+ endif
+#ifdef MPI
+ if (Kwita.eq.0) call recv_stop_sig(Kwita)
+ if (Kwita.lt.0) then
+ write (iout,*) 'Stop signal received. Terminating.'
+ write (*,*) 'Stop signal received. Terminating.'
+ nthread=ithread
+ write (*,*) 'nthread=',nthread,' ithread=',ithread
+ goto 777
+ endif
+#endif
+ enddo
+#ifdef MPI
+ call send_stop_sig(-1)
+#endif
+ 777 continue
+#ifdef MPI
+C Any messages left for me?
+ call pattern_receive
+ if (Kwita.eq.0) call recv_stop_sig(Kwita)
+#endif
+ call write_thread_summary
+#ifdef MPI
+ if (king.eq.king) then
+ Kwita=1
+ do while (Kwita.ne.0 .or. nacc_tot.ne.0)
+ Kwita=0
+ nacc_tot=0
+ call recv_stop_sig(Kwita)
+ call receive_MCM_info
+ enddo
+ do iproc=1,nprocs-1
+ call receive_thread_results(iproc)
+ enddo
+ call write_thread_summary
+ else
+ call send_thread_results
+ endif
+#endif
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine write_thread_summary
+C Thread the sequence through a database of known structures
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DBASE'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.THREAD'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.HEADER'
+ include 'COMMON.NAMES'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.TIME1'
+#ifdef MPI
+ include 'COMMON.INFO'
+#endif
+ dimension ip(maxthread)
+ double precision energia(0:n_ene)
+ write (iout,'(30x,a/)')
+ & ' *********** Summary threading statistics ************'
+ write (iout,'(a)') 'Initial energies:'
+ write (iout,'(a4,2x,a12,14a14,3a8)')
+ & 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',
+ & 'RMSnat','NatCONT','NNCONT','RMS'
+C Energy sort patterns
+ do i=1,nthread
+ ip(i)=i
+ enddo
+ do i=1,nthread-1
+ enet=ener(n_ene-1,ip(i))
+ jj=i
+ do j=i+1,nthread
+ if (ener(n_ene-1,ip(j)).lt.enet) then
+ jj=j
+ enet=ener(n_ene-1,ip(j))
+ endif
+ enddo
+ if (jj.ne.i) then
+ ipj=ip(jj)
+ ip(jj)=ip(i)
+ ip(i)=ipj
+ endif
+ enddo
+ do ik=1,nthread
+ i=ip(ik)
+ ii=ipatt(1,i)
+ ist=nres_base(2,ii)+ipatt(2,i)
+ do kk=1,n_ene_comp
+ energia(i)=ener0(kk,i)
+ enddo
+ etot=ener0(n_ene_comp+1,i)
+ rmsnat=ener0(n_ene_comp+2,i)
+ rms=ener0(n_ene_comp+3,i)
+ frac=ener0(n_ene_comp+4,i)
+ frac_nn=ener0(n_ene_comp+5,i)
+
+ if (refstr) then
+ write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)')
+ & i,str_nam(ii),ist+1,
+ & (energia(print_order(kk)),kk=1,nprint_ene),
+ & etot,rmsnat,frac,frac_nn,rms
+ else
+ write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)')
+ & i,str_nam(ii),ist+1,
+ & (energia(print_order(kk)),kk=1,nprint_ene),etot
+ endif
+ enddo
+ write (iout,'(//a)') 'Final energies:'
+ write (iout,'(a4,2x,a12,17a14,3a8)')
+ & 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',
+ & 'RMSnat','NatCONT','NNCONT','RMS'
+ do ik=1,nthread
+ i=ip(ik)
+ ii=ipatt(1,i)
+ ist=nres_base(2,ii)+ipatt(2,i)
+ do kk=1,n_ene_comp
+ energia(kk)=ener(kk,ik)
+ enddo
+ etot=ener(n_ene_comp+1,i)
+ rmsnat=ener(n_ene_comp+2,i)
+ rms=ener(n_ene_comp+3,i)
+ frac=ener(n_ene_comp+4,i)
+ frac_nn=ener(n_ene_comp+5,i)
+ write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)')
+ & i,str_nam(ii),ist+1,
+ & (energia(print_order(kk)),kk=1,nprint_ene),
+ & etot,rmsnat,frac,frac_nn,rms
+ enddo
+ write (iout,'(/a/)') 'IEXAM array:'
+ write (iout,'(i5)') nexcl
+ do i=1,nexcl
+ write (iout,'(2i5)') iexam(1,i),iexam(2,i)
+ enddo
+ write (iout,'(/a,1pe14.4/a,1pe14.4/)')
+ & 'Max. time for threading step ',max_time_for_thread,
+ & 'Average time for threading step: ',ave_time_for_thread
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine sc_conf
+C Sample (hopefully) optimal SC orientations given backcone conformation.
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DBASE'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.THREAD'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.HEADER'
+ include 'COMMON.GEO'
+ include 'COMMON.IOUNITS'
+ double precision varia(maxvar)
+ common /srutu/ icall
+ double precision energia(0:n_ene)
+ logical glycine,fail
+ maxsample=10
+ link_end0=link_end
+ link_end=min0(link_end,nss)
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+cd print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1)
+ call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail)
+ endif
+ enddo
+ call chainbuild
+ call etotal(energia(0))
+ do isample=1,maxsample
+C Choose a non-glycine side chain.
+ glycine=.true.
+ do while(glycine)
+ ind_sc=iran_num(nnt,nct)
+ glycine=(itype(ind_sc).eq.10)
+ enddo
+ alph0=alph(ind_sc)
+ omeg0=omeg(ind_sc)
+ call gen_side(itype(ind_sc),theta(ind_sc+1),alph(ind_sc),
+ & omeg(ind_sc),fail)
+ call chainbuild
+ call etotal(energia(0))
+cd write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))')
+cd & 'Step:',isample,' SC',ind_sc,' alpha',alph(ind_sc)*rad2deg,
+cd & ' omega',omeg(ind_sc)*rad2deg,' old energy',e0,' new energy',e1
+ e1=0.0d0
+ if (e0.le.e1) then
+ alph(ind_sc)=alph0
+ omeg(ind_sc)=omeg0
+ else
+ e0=e1
+ endif
+ enddo
+ link_end=link_end0
+ return
+ end
+c---------------------------------------------------------------------------
+ subroutine write_stat_thread(ithread,ipattern,ist)
+ implicit real*8 (a-h,o-z)
+ include "DIMENSIONS"
+ include "COMMON.CONTROL"
+ include "COMMON.IOUNITS"
+ include "COMMON.THREAD"
+ include "COMMON.FFIELD"
+ include "COMMON.DBASE"
+ include "COMMON.NAMES"
+ double precision energia(0:n_ene)
+
+#if defined(AIX) || defined(PGI)
+ open(istat,file=statname,position='append')
+#else
+ open(istat,file=statname,access='append')
+#endif
+ do i=1,n_ene_comp
+ energia(i)=ener(i,ithread)
+ enddo
+ etot=ener(n_ene_comp+1,ithread)
+ rmsnat=ener(n_ene_comp+2,ithread)
+ rms=ener(n_ene_comp+3,ithread)
+ frac=ener(n_ene_comp+4,ithread)
+ frac_nn=ener(n_ene_comp+5,ithread)
+ write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)')
+ & ithread,str_nam(ipattern),ist+1,
+ & (energia(print_order(i)),i=1,nprint_ene),
+ & etot,rmsnat,frac,frac_nn,rms
+ close (istat)
+ return
+ end
projects / unres.git / blobdiff