Merge branch 'bartek' of mmka.chem.univ.gda.pl:unres into bartek
[unres.git] / 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