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