ctest wham & cluster src-M
[unres.git] / source / unres / src_CSA / unres_csa.F
1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2 C                                                                              C
3 C                                U N R E S                                     C
4 C                                                                              C
5 C Program to carry out conformational search of proteins in an united-residue  C
6 C approximation.                                                               C
7 C                                                                              C
8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
9       implicit real*8 (a-h,o-z)
10       include 'DIMENSIONS'
11
12
13 #ifdef MPI
14       include 'mpif.h'
15       include 'COMMON.SETUP'
16 #endif
17       include 'COMMON.TIME1'
18       include 'COMMON.INTERACT'
19       include 'COMMON.NAMES'
20       include 'COMMON.GEO'
21       include 'COMMON.HEADER'
22       include 'COMMON.CONTROL'
23       include 'COMMON.CONTACTS'
24       include 'COMMON.CHAIN'
25       include 'COMMON.VAR'
26       include 'COMMON.IOUNITS'
27       include 'COMMON.FFIELD'
28 c      include 'COMMON.REMD'
29 c      include 'COMMON.MD'
30       include 'COMMON.SBRIDGE'
31       double precision hrtime,mintime,sectime
32       character*64 text_mode_calc(-2:14) /'test',
33      & 'SC rotamer distribution',
34      & 'Energy evaluation or minimization',
35      & 'Regularization of PDB structure',
36      & 'Threading of a sequence on PDB structures',
37      & 'Monte Carlo (with minimization) ',
38      & 'Energy minimization of multiple conformations',
39      & 'Checking energy gradient',
40      & 'Entropic sampling Monte Carlo (with minimization)',
41      & 'Energy map',
42      & 'CSA calculations',
43      & 'Not used 9',
44      & 'Not used 10',
45      & 'Soft regularization of PDB structure',
46      & 'Mesoscopic molecular dynamics (MD) ',
47      & 'Not used 13',
48      & 'Replica exchange molecular dynamics (REMD)'/
49       external ilen
50
51 c      call memmon_print_usage()
52
53       call init_task
54       if (me.eq.king)
55      & write(iout,*)'### LAST MODIFIED  11/03/09 1:19PM by czarek'  
56       if (me.eq.king) call cinfo
57 C Read force field parameters and job setup data
58       call readrtns
59       call flush(iout)
60 C
61       if (me.eq.king .or. .not. out1file) then
62        write (iout,'(2a/)') 
63      & text_mode_calc(modecalc)(:ilen(text_mode_calc(modecalc))),
64      & ' calculation.' 
65        if (minim) write (iout,'(a)') 
66      &  'Conformations will be energy-minimized.'
67        write (iout,'(80(1h*)/)') 
68       endif
69       call flush(iout)
70 C
71 c      if (modecalc.eq.-2) then
72 c        call test
73 c        stop
74 c      else if (modecalc.eq.-1) then
75 c        write(iout,*) "call check_sc_map next"
76 c        call check_bond
77 c        stop
78 c      endif
79 #ifdef MPI
80       if (fg_rank.gt.0) then
81 C Fine-grain slaves just do energy and gradient components.
82         call ergastulum ! slave workhouse in Latin
83       else
84 #endif
85       if (modecalc.eq.0) then
86         call exec_eeval_or_minim
87 c      else if (modecalc.eq.1) then
88 c        call exec_regularize
89 c      else if (modecalc.eq.2) then
90 c        call exec_thread
91 c      else if (modecalc.eq.3 .or. modecalc .eq.6) then
92 c        call exec_MC
93       else if (modecalc.eq.4) then
94         call exec_mult_eeval_or_minim
95       else if (modecalc.eq.5) then
96          call exec_checkgrad
97 c      else if (ModeCalc.eq.7) then
98 c        call exec_map
99       else if (ModeCalc.eq.8) then
100         call exec_CSA
101 c      else if (modecalc.eq.11) then
102 c        call exec_softreg
103 c      else if (modecalc.eq.12) then
104 c        call exec_MD
105 c      else if (modecalc.eq.14) then
106 c        call exec_MREMD
107       else
108         write (iout,'(a)') 'This calculation type is not supported',
109      &   ModeCalc
110       endif
111 #ifdef MPI
112       endif
113 C Finish task.
114       if (fg_rank.eq.0) call finish_task
115 c      call memmon_print_usage()
116 #ifdef TIMING
117        call print_detailed_timing
118 #endif
119       call MPI_Finalize(ierr)
120       stop 'Bye Bye...'
121 #else
122       call dajczas(tcpu(),hrtime,mintime,sectime)
123       stop '********** Program terminated normally.'
124 #endif
125       end
126 c---------------------------------------------------------------------------
127       subroutine exec_eeval_or_minim
128       implicit real*8 (a-h,o-z)
129       include 'DIMENSIONS'
130 #ifdef MPI
131       include 'mpif.h'
132 #endif
133       include 'COMMON.SETUP'
134       include 'COMMON.TIME1'
135       include 'COMMON.INTERACT'
136       include 'COMMON.NAMES'
137       include 'COMMON.GEO'
138       include 'COMMON.HEADER'
139       include 'COMMON.CONTROL'
140       include 'COMMON.CONTACTS'
141       include 'COMMON.CHAIN'
142       include 'COMMON.VAR'
143       include 'COMMON.IOUNITS'
144       include 'COMMON.FFIELD'
145 c      include 'COMMON.REMD'
146 c      include 'COMMON.MD'
147       include 'COMMON.SBRIDGE'
148       common /srutu/ icall
149       double precision energy(0:n_ene)
150       double precision energy_long(0:n_ene),energy_short(0:n_ene)
151       double precision varia(maxvar)
152       if (indpdb.eq.0) call chainbuild
153 c      time00=MPI_Wtime()
154       call chainbuild_cart
155 c      if (split_ene) then
156 c       print *,"Processor",myrank," after chainbuild"
157 c       icall=1
158 c       call etotal_long(energy_long(0))
159 c       write (iout,*) "Printing long range energy"
160 c       call enerprint(energy_long(0))
161 c       call etotal_short(energy_short(0))
162 c       write (iout,*) "Printing short range energy"
163 c       call enerprint(energy_short(0))
164 c       do i=0,n_ene
165 c         energy(i)=energy_long(i)+energy_short(i)
166 c         write (iout,*) i,energy_long(i),energy_short(i),energy(i)
167 c       enddo
168 c       write (iout,*) "Printing long+short range energy"
169 c       call enerprint(energy(0))
170 c      endif
171       call etotal(energy(0))
172 c      time_ene=MPI_Wtime()-time00
173       write (iout,*) "Time for energy evaluation",time_ene
174       print *,"after etotal"
175       etota = energy(0)
176       etot =etota
177       call enerprint(energy(0))
178 c      call hairpin(.true.,nharp,iharp)
179 c      call secondary2(.true.)
180       if (minim) then
181
182         if (dccart) then
183           print *, 'Calling MINIM_DC'
184 c          time1=MPI_WTIME()
185           call minim_dc(etot,iretcode,nfun)
186         else
187           if (indpdb.ne.0) then 
188             call bond_regular
189             call chainbuild
190           endif
191           call geom_to_var(nvar,varia)
192           print *,'Calling MINIMIZE.'
193 c          time1=MPI_WTIME()
194           call minimize(etot,varia,iretcode,nfun)
195         endif
196         print *,'SUMSL return code is',iretcode,' eval ',nfun
197 c        evals=nfun/(MPI_WTIME()-time1)
198         print *,'# eval/s',evals
199         print *,'refstr=',refstr
200 c        call hairpin(.true.,nharp,iharp)
201 c        call secondary2(.true.)
202         call etotal(energy(0))
203         etot = energy(0)
204         call enerprint(energy(0))
205
206         call intout
207         call briefout(0,etot)
208         if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
209           write (iout,'(a,i3)') 'SUMSL return code:',iretcode
210           write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
211           write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
212       else
213         print *,'refstr=',refstr
214         if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
215         call briefout(0,etot)
216       endif
217       if (outpdb) call pdbout(etot,titel(:32),ipdb)
218       if (outmol2) call mol2out(etot,titel(:32))
219       return
220       end
221
222       subroutine exec_checkgrad
223       implicit real*8 (a-h,o-z)
224       include 'DIMENSIONS'
225 #ifdef MPI
226       include 'mpif.h'
227 #endif
228       include 'COMMON.SETUP'
229       include 'COMMON.TIME1'
230       include 'COMMON.INTERACT'
231       include 'COMMON.NAMES'
232       include 'COMMON.GEO'
233       include 'COMMON.HEADER'
234       include 'COMMON.CONTROL'
235       include 'COMMON.CONTACTS'
236       include 'COMMON.CHAIN'
237       include 'COMMON.VAR'
238       include 'COMMON.IOUNITS'
239       include 'COMMON.FFIELD'
240 c      include 'COMMON.REMD'
241       include 'COMMON.MD_'
242       include 'COMMON.SBRIDGE'
243       common /srutu/ icall
244       double precision energy(0:max_ene)
245 c      do i=2,nres
246 c        vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0)
247 c        if (itype(i).ne.10) 
248 c     &      vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0)
249 c      enddo
250       if (indpdb.eq.0) call chainbuild
251 c      do i=0,nres
252 c        do j=1,3
253 c          dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0)
254 c        enddo
255 c      enddo
256 c      do i=1,nres-1
257 c        if (itype(i).ne.10) then
258 c          do j=1,3
259 c            dc(j,i+nres)=dc(j,i+nres)+ran_number(-0.2d0,0.2d0)
260 c          enddo
261 c        endif
262 c      enddo
263 c      do j=1,3
264 c        dc(j,0)=ran_number(-0.2d0,0.2d0)
265 c      enddo
266       usampl=.true.
267       totT=1.d0
268       eq_time=0.0d0
269 c      call read_fragments
270 c      read(inp,*) t_bath
271 c      call rescale_weights(t_bath)
272       call chainbuild_cart
273       call cartprint
274       call intout
275       icall=1
276       call etotal(energy(0))
277       etot = energy(0)
278       call enerprint(energy(0))
279       write (iout,*) "Uconst",Uconst," Uconst_back",uconst_back
280       print *,'icheckgrad=',icheckgrad
281       goto (10,20,30) icheckgrad
282   10  call check_ecartint
283       return
284   20  call check_cartgrad
285       return
286   30  call check_eint
287       return
288       end
289 c---------------------------------------------------------------------------
290       subroutine exec_CSA
291 #ifdef MPI
292       include "mpif.h"
293 #endif
294       include 'DIMENSIONS'
295       include 'COMMON.IOUNITS'
296 C Conformational Space Annealling programmed by Jooyoung Lee.
297 C This method works only with parallel machines!
298 #ifdef MPI
299       call together
300 #else
301       write (iout,*) "CSA works on parallel machines only"
302 #endif
303       return
304       end
305 c---------------------------------------------------------------------------
306 #ifdef MPI
307       subroutine exec_mult_eeval_or_minim
308       implicit real*8 (a-h,o-z)
309       include 'DIMENSIONS'
310       include 'mpif.h'
311       integer muster(mpi_status_size)
312       include 'COMMON.SETUP'
313       include 'COMMON.TIME1'
314       include 'COMMON.INTERACT'
315       include 'COMMON.NAMES'
316       include 'COMMON.GEO'
317       include 'COMMON.HEADER'
318       include 'COMMON.CONTROL'
319       include 'COMMON.CONTACTS'
320       include 'COMMON.CHAIN'
321       include 'COMMON.VAR'
322       include 'COMMON.IOUNITS'
323       include 'COMMON.FFIELD'
324       include 'COMMON.SBRIDGE'
325       double precision varia(maxvar)
326       integer ind(6)
327       double precision energy(0:n_ene)
328       logical eof
329       eof=.false.
330
331       if(me.ne.king) then
332         call minim_mcmf
333         return
334       endif
335
336       close (intin)
337       open(intin,file=intinname,status='old')
338       write (istat,'(a5,100a12)')"#    ",
339      &  (wname(print_order(i)),i=1,nprint_ene)
340       if (refstr) then
341         write (istat,'(a5,100a12)')"#    ",
342      &   (ename(print_order(i)),i=1,nprint_ene),
343      &   "ETOT total","RMSD","nat.contact","nnt.contact",
344      &   "cont.order","TMscore"
345       else
346         write (istat,'(a5,100a12)')"#    ",
347      &    (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
348       endif
349
350       if (.not.minim) then
351         do while (.not. eof)
352           if (read_cart) then
353             read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene
354             call read_x(intin,*11)
355 c Broadcast the order to compute internal coordinates to the slaves.
356             if (nfgtasks.gt.1)
357      &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
358             call int_from_cart1(.false.)
359           else
360             read (intin,'(i5)',end=1100,err=1100) iconf
361             call read_angles(intin,*11)
362             call geom_to_var(nvar,varia)
363             call chainbuild
364           endif
365           write (iout,'(a,i7)') 'Conformation #',iconf
366           call etotal(energy(0))
367           call briefout(iconf,energy(0))
368           call enerprint(energy(0))
369           etot=energy(0)
370           if (refstr) then 
371             call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
372             call calc_tmscore(tm,.true.)
373             write (istat,'(i5,100(f12.3))') iconf,
374      &      (energy(print_order(i)),i=1,nprint_ene),etot,
375      &       rms,frac,frac_nn,co,tm
376           else
377             write (istat,'(i5,100(f12.3))') iconf,
378      &     (energy(print_order(i)),i=1,nprint_ene),etot
379           endif
380         enddo
381 1100    continue
382         goto 1101
383       endif
384
385       mm=0
386       imm=0
387       nft=0
388       ene0=0.0d0
389       n=0
390       iconf=0
391       do while (.not. eof)
392         mm=mm+1
393         if (mm.lt.nodes) then
394           if (read_cart) then
395             read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
396             call read_x(intin,*11)
397 c Broadcast the order to compute internal coordinates to the slaves.
398             if (nfgtasks.gt.1) 
399      &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
400             call int_from_cart1(.false.)
401           else
402             read (intin,'(i5)',end=11,err=11) iconf
403             call read_angles(intin,*11)
404             call geom_to_var(nvar,varia)
405             call chainbuild
406           endif
407
408           n=n+1
409           write (iout,*) 'Conformation #',iconf,' read'
410          imm=imm+1
411          ind(1)=1
412          ind(2)=n
413          ind(3)=0
414          ind(4)=0
415          ind(5)=0
416          ind(6)=0
417          ene0=0.0d0
418          call mpi_send(ind,6,mpi_integer,mm,idint,CG_COMM,
419      *                  ierr)
420          call mpi_send(varia,nvar,mpi_double_precision,mm,
421      *                  idreal,CG_COMM,ierr)
422          call mpi_send(ene0,1,mpi_double_precision,mm,
423      *                  idreal,CG_COMM,ierr)
424 c         print *,'task ',n,' sent to worker ',mm,nvar
425         else
426          call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
427      *                 CG_COMM,muster,ierr)
428          man=muster(mpi_source)
429 c         print *,'receiving result from worker ',man,' (',iii1,iii,')'
430          call mpi_recv(varia,nvar,mpi_double_precision, 
431      *               man,idreal,CG_COMM,muster,ierr)
432          call mpi_recv(ene,1,
433      *               mpi_double_precision,man,idreal,
434      *               CG_COMM,muster,ierr)
435          call mpi_recv(ene0,1,
436      *               mpi_double_precision,man,idreal,
437      *               CG_COMM,muster,ierr)
438 c         print *,'result received from worker ',man,' sending now'
439
440           call var_to_geom(nvar,varia)
441           call chainbuild
442           call etotal(energy(0))
443           iconf=ind(2)
444           write (iout,*)
445           write (iout,*)
446           write (iout,*) 'Conformation #',iconf," sumsl return code ",
447      &                      ind(5)
448
449           etot=energy(0)
450           call enerprint(energy(0))
451           call briefout(iconf,etot)
452           if (refstr) then 
453             call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
454             call calc_tmscore(tm,.true.)
455             write (istat,'(i5,100(f12.3))') iconf,
456      &     (energy(print_order(i)),i=1,nprint_ene),etot,
457      &     rms,frac,frac_nn,co,tm
458           else
459             write (istat,'(i5,100(f12.3))') iconf,
460      &     (energy(print_order(i)),i=1,nprint_ene),etot
461           endif
462
463           imm=imm-1
464           if (read_cart) then
465             read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
466             call read_x(intin,*11)
467 c Broadcast the order to compute internal coordinates to the slaves.
468             if (nfgtasks.gt.1)
469      &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
470             call int_from_cart1(.false.)
471           else
472             read (intin,'(i5)',end=11,err=11) iconf
473             call read_angles(intin,*11)
474             call geom_to_var(nvar,varia)
475             call chainbuild
476           endif
477           n=n+1
478           write (iout,*) 'Conformation #',iconf,' read'
479           imm=imm+1
480           ind(1)=1
481           ind(2)=n
482           ind(3)=0
483           ind(4)=0
484           ind(5)=0
485           ind(6)=0
486           call mpi_send(ind,6,mpi_integer,man,idint,CG_COMM,
487      *                  ierr)
488           call mpi_send(varia,nvar,mpi_double_precision,man, 
489      *                  idreal,CG_COMM,ierr)
490           call mpi_send(ene0,1,mpi_double_precision,man,
491      *                  idreal,CG_COMM,ierr)
492           nf_mcmf=nf_mcmf+ind(4)
493           nmin=nmin+1
494         endif
495       enddo
496 11    continue
497       do j=1,imm
498         call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
499      *               CG_COMM,muster,ierr)
500         man=muster(mpi_source)
501         call mpi_recv(varia,nvar,mpi_double_precision, 
502      *               man,idreal,CG_COMM,muster,ierr)
503         call mpi_recv(ene,1,
504      *               mpi_double_precision,man,idreal,
505      *               CG_COMM,muster,ierr)
506         call mpi_recv(ene0,1,
507      *               mpi_double_precision,man,idreal,
508      *               CG_COMM,muster,ierr)
509
510         call var_to_geom(nvar,varia)
511         call chainbuild
512         call etotal(energy(0))
513         iconf=ind(2)
514         write (iout,*)
515         write (iout,*)
516         write (iout,*) 'Conformation #',iconf," sumsl return code ",
517      &                  ind(5)
518
519         etot=energy(0)
520         call enerprint(energy(0))
521         call briefout(iconf,etot)
522         if (refstr) then 
523           call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
524           call calc_tmscore(tm,.true.)
525           write (istat,'(i5,100(f12.3))') iconf,
526      &   (energy(print_order(i)),i=1,nprint_ene),etot,
527      &   rms,frac,frac_nn,co,tm
528         else
529           write (istat,'(i5,100(f12.3))') iconf,
530      &    (energy(print_order(i)),i=1,nprint_ene),etot
531         endif
532         nmin=nmin+1
533       enddo
534 1101  continue
535       do i=1, nodes-1
536          ind(1)=0
537          ind(2)=0
538          ind(3)=0
539          ind(4)=0
540          ind(5)=0
541          ind(6)=0
542          call mpi_send(ind,6,mpi_integer,i,idint,CG_COMM,
543      *                  ierr)
544       enddo
545       return
546       end
547 #else
548       subroutine exec_mult_eeval_or_minim
549       include 'DIMENSIONS'
550       include 'COMMON.IOUNITS'
551       write (iout,*) "Unsupported option in serial version"
552       return
553       end
554 #endif
555 c---------------------------------------------------------------------------
556