X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fthread.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fthread.F;h=f7137443f82ce7814b2230434d50531e479b55ed;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/thread.F b/source/unres/src_MD-M-newcorr/thread.F new file mode 100644 index 0000000..f713744 --- /dev/null +++ b/source/unres/src_MD-M-newcorr/thread.F @@ -0,0 +1,549 @@ + 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.ntyp1) 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.ntyp1) 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