ctest wham & cluster src-M
[unres.git] / source / unres / src_MD / minimize_p.F
1       subroutine minimize(etot,x,iretcode,nfun)
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
5 *********************************************************************
6 * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
7 * the calling subprogram.                                           *     
8 * when d(i)=1.0, then v(35) is the length of the initial step,      *     
9 * calculated in the usual pythagorean way.                          *     
10 * absolute convergence occurs when the function is within v(31) of  *     
11 * zero. unless you know the minimum value in advance, abs convg     *     
12 * is probably not useful.                                           *     
13 * relative convergence is when the model predicts that the function *   
14 * will decrease by less than v(32)*abs(fun).                        *   
15 *********************************************************************
16       include 'COMMON.IOUNITS'
17       include 'COMMON.VAR'
18       include 'COMMON.GEO'
19       include 'COMMON.MINIM'
20       common /srutu/ icall
21       dimension iv(liv)                                               
22       double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
23       double precision energia(0:n_ene)
24       external func,gradient,fdum
25       external func_restr,grad_restr
26       logical not_done,change,reduce 
27 c      common /przechowalnia/ v
28
29       icall = 1
30
31       NOT_DONE=.TRUE.
32
33 c     DO WHILE (NOT_DONE)
34
35       call deflt(2,iv,liv,lv,v)                                         
36 * 12 means fresh start, dont call deflt                                 
37       iv(1)=12                                                          
38 * max num of fun calls                                                  
39       if (maxfun.eq.0) maxfun=500
40       iv(17)=maxfun
41 * max num of iterations                                                 
42       if (maxmin.eq.0) maxmin=1000
43       iv(18)=maxmin
44 * controls output                                                       
45       iv(19)=2                                                          
46 * selects output unit                                                   
47       iv(21)=0
48       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
49 * 1 means to print out result                                           
50       iv(22)=print_min_res
51 * 1 means to print out summary stats                                    
52       iv(23)=print_min_stat
53 * 1 means to print initial x and d                                      
54       iv(24)=print_min_ini
55 * min val for v(radfac) default is 0.1                                  
56       v(24)=0.1D0                                                       
57 * max val for v(radfac) default is 4.0                                  
58       v(25)=2.0D0                                                       
59 c     v(25)=4.0D0                                                       
60 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
61 * the sumsl default is 0.1                                              
62       v(26)=0.1D0
63 * false conv if (act fnctn decrease) .lt. v(34)                         
64 * the sumsl default is 100*machep                                       
65       v(34)=v(34)/100.0D0                                               
66 * absolute convergence                                                  
67       if (tolf.eq.0.0D0) tolf=1.0D-4
68       v(31)=tolf
69 * relative convergence                                                  
70       if (rtolf.eq.0.0D0) rtolf=1.0D-4
71       v(32)=rtolf
72 * controls initial step size                                            
73        v(35)=1.0D-1                                                    
74 * large vals of d correspond to small components of step                
75       do i=1,nphi
76         d(i)=1.0D-1
77       enddo
78       do i=nphi+1,nvar
79         d(i)=1.0D-1
80       enddo
81 cd    print *,'Calling SUMSL'
82 c     call var_to_geom(nvar,x)
83 c     call chainbuild
84 c     call etotal(energia(0))
85 c     etot = energia(0)
86       IF (mask_r) THEN
87        call x2xx(x,xx,nvar_restr)
88        call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
89      &                    iv,liv,lv,v,idum,rdum,fdum)      
90        call xx2x(x,xx)
91       ELSE
92        call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)      
93       ENDIF
94       etot=v(10)                                                      
95       iretcode=iv(1)
96 cd    print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
97 cd    write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
98 c     call intout
99 c     change=reduce(x)
100       call var_to_geom(nvar,x)
101 c     if (change) then
102 c       write (iout,'(a)') 'Reduction worked, minimizing again...'
103 c     else
104 c       not_done=.false.
105 c     endif
106       call chainbuild
107 c     call etotal(energia(0))
108 c     etot=energia(0)
109 c     call enerprint(energia(0))
110       nfun=iv(6)
111
112 c     write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
113
114 c     ENDDO ! NOT_DONE
115
116       return  
117       end  
118 #ifdef MPI
119 c----------------------------------------------------------------------------
120       subroutine ergastulum
121       implicit real*8 (a-h,o-z)
122       include 'DIMENSIONS'
123 #ifdef MPI
124       include "mpif.h"
125 #endif
126       include 'COMMON.SETUP'
127       include 'COMMON.DERIV'
128       include 'COMMON.VAR'
129       include 'COMMON.IOUNITS'
130       include 'COMMON.FFIELD'
131       include 'COMMON.INTERACT'
132       include 'COMMON.MD'
133       include 'COMMON.TIME1'
134       double precision z(maxres6),d_a_tmp(maxres6)
135       double precision edum(0:n_ene),time_order(0:10)
136       double precision Gcopy(maxres2,maxres2)
137       double precision e_tmp
138       common /przechowalnia/ Gcopy
139       integer icall /0/
140 C Workers wait for variables and NF, and NFL from the boss 
141       iorder=0
142       do while (iorder.ge.0)
143 c      write (*,*) 'Processor',fg_rank,' CG group',kolor,
144 c     & ' receives order from Master'
145         time00=MPI_Wtime()
146         call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
147         time_Bcast=time_Bcast+MPI_Wtime()-time00
148         if (icall.gt.4 .and. iorder.ge.0) 
149      &   time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
150        icall=icall+1
151 c      write (*,*) 
152 c     & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
153         if (iorder.eq.0) then
154           call zerograd
155           call etotal(edum)
156 c          write (2,*) "After etotal"
157 c          write (2,*) "dimen",dimen," dimen3",dimen3
158 c          call flush(2)
159         else if (iorder.eq.2) then
160           call zerograd
161           call etotal_short(edum)
162 c          write (2,*) "After etotal_short"
163 c          write (2,*) "dimen",dimen," dimen3",dimen3
164 c          call flush(2)
165         else if (iorder.eq.3) then
166           call zerograd
167           call etotal_long(edum)
168 c          write (2,*) "After etotal_long"
169 c          write (2,*) "dimen",dimen," dimen3",dimen3
170 c          call flush(2)
171         else if (iorder.eq.1) then
172           call sum_gradient
173 c          write (2,*) "After sum_gradient"
174 c          write (2,*) "dimen",dimen," dimen3",dimen3
175 c          call flush(2)
176         else if (iorder.eq.4) then
177           call ginv_mult(z,d_a_tmp)
178         else if (iorder.eq.5) then
179 c Setup MD things for a slave
180           dimen=(nct-nnt+1)+nside
181           dimen1=(nct-nnt)+(nct-nnt+1)
182           dimen3=dimen*3
183 c          write (2,*) "dimen",dimen," dimen3",dimen3
184 c          call flush(2)
185           call int_bounds(dimen,igmult_start,igmult_end)
186           igmult_start=igmult_start-1
187           call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
188      &     ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
189            my_ng_count=igmult_end-igmult_start
190           call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
191      &     MPI_INTEGER,FG_COMM,IERROR)
192 c          write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
193 c          write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
194           myginv_ng_count=maxres2*my_ng_count
195 c          write (2,*) "igmult_start",igmult_start," igmult_end",
196 c     &     igmult_end," my_ng_count",my_ng_count
197 c          call flush(2)
198           call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
199      &     nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
200           call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
201      &     nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
202 c          write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
203 c          write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
204 c          call flush(2)
205 c          call MPI_Barrier(FG_COMM,IERROR)
206           time00=MPI_Wtime()
207           call MPI_Scatterv(ginv(1,1),nginv_counts(0),
208      &     nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
209      &     myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
210 #ifdef TIMING
211           time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
212 #endif
213           do i=1,dimen
214             do j=1,2*my_ng_count
215               ginv(j,i)=gcopy(i,j)
216             enddo
217           enddo
218 c          write (2,*) "dimen",dimen," dimen3",dimen3
219 c          write (2,*) "End MD setup"
220 c          call flush(2)
221 c           write (iout,*) "My chunk of ginv_block"
222 c           call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
223         else if (iorder.eq.6) then
224           call int_from_cart1(.false.)
225         else if (iorder.eq.7) then
226           call chainbuild_cart
227         else if (iorder.eq.8) then
228           call intcartderiv
229         else if (iorder.eq.9) then
230           call fricmat_mult(z,d_a_tmp)
231         else if (iorder.eq.10) then
232           call setup_fricmat
233         else if (iorder.eq.11) then
234 c broadcast new value of iset
235           call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
236         else if (iorder.eq.12) then
237 c broadcast new value of iset and
238 c calculate restraint homology energy and
239 c sum it over FG_COMM
240           call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
241           call e_modeller(e_tmp)
242           call MPI_Barrier(FG_COMM,IERR)
243           call MPI_Reduce(e_tmp,ehomology_constr,1,
244      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
245         else if (iorder.eq.13) then
246           call chainbuild
247         endif
248       enddo
249       write (*,*) 'Processor',fg_rank,' CG group',kolor,
250      &  ' absolute rank',myrank,' leves ERGASTULUM.'
251       write(*,*)'Processor',fg_rank,' wait times for respective orders',
252      & (' order[',i,']',time_order(i),i=0,12)
253       return
254       end
255 #endif
256 ************************************************************************
257       subroutine func(n,x,nf,f,uiparm,urparm,ufparm)  
258       implicit real*8 (a-h,o-z)
259       include 'DIMENSIONS'
260       include 'COMMON.DERIV'
261       include 'COMMON.IOUNITS'
262       include 'COMMON.GEO'
263       common /chuju/ jjj
264       double precision energia(0:n_ene)
265       integer jjj
266       double precision ufparm
267       external ufparm                                                   
268       integer uiparm(1)                                                 
269       real*8 urparm(1)                                                    
270       dimension x(maxvar)
271 c     if (jjj.gt.0) then
272 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
273 c     endif
274       nfl=nf
275       icg=mod(nf,2)+1
276 cd      print *,'func',nf,nfl,icg
277       call var_to_geom(n,x)
278       call zerograd
279       call chainbuild
280 cd    write (iout,*) 'ETOTAL called from FUNC'
281       call etotal(energia(0))
282       call sum_gradient
283       f=energia(0)
284 c     if (jjj.gt.0) then
285 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
286 c       write (iout,*) 'f=',f
287 c       jjj=0
288 c     endif
289       return                                                            
290       end                                                               
291 ************************************************************************
292       subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)  
293       implicit real*8 (a-h,o-z)
294       include 'DIMENSIONS'
295       include 'COMMON.DERIV'
296       include 'COMMON.IOUNITS'
297       include 'COMMON.GEO'
298       common /chuju/ jjj
299       double precision energia(0:n_ene)
300       integer jjj
301       double precision ufparm
302       external ufparm                                                   
303       integer uiparm(1)                                                 
304       real*8 urparm(1)                                                    
305       dimension x(maxvar)
306 c     if (jjj.gt.0) then
307 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
308 c     endif
309       nfl=nf
310       icg=mod(nf,2)+1
311       call var_to_geom_restr(n,x)
312       call zerograd
313       call chainbuild
314 cd    write (iout,*) 'ETOTAL called from FUNC'
315       call etotal(energia(0))
316       call sum_gradient
317       f=energia(0)
318 c     if (jjj.gt.0) then
319 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
320 c       write (iout,*) 'f=',etot
321 c       jjj=0
322 c     endif
323       return                                                            
324       end                                                               
325 c-------------------------------------------------------
326       subroutine x2xx(x,xx,n)
327       implicit real*8 (a-h,o-z)
328       include 'DIMENSIONS'
329       include 'COMMON.VAR'
330       include 'COMMON.CHAIN'
331       include 'COMMON.INTERACT'
332       double precision xx(maxvar),x(maxvar)
333
334       do i=1,nvar
335         varall(i)=x(i)
336       enddo
337
338       ig=0                                                                      
339       igall=0                                                                   
340       do i=4,nres                                                               
341         igall=igall+1                                                           
342         if (mask_phi(i).eq.1) then                                              
343           ig=ig+1                                                               
344           xx(ig)=x(igall)                       
345         endif                                                                   
346       enddo                                                                     
347                                                                                 
348       do i=3,nres                                                               
349         igall=igall+1                                                           
350         if (mask_theta(i).eq.1) then                                            
351           ig=ig+1                                                               
352           xx(ig)=x(igall)
353         endif                                                                   
354       enddo                                          
355
356       do ij=1,2                                                                 
357       do i=2,nres-1                                                             
358         if (itype(i).ne.10) then                                                
359           igall=igall+1                                                         
360           if (mask_side(i).eq.1) then                                           
361             ig=ig+1                                                             
362             xx(ig)=x(igall)
363           endif                                                                 
364         endif                                                                   
365       enddo                                                                     
366       enddo                              
367  
368       n=ig
369
370       return
371       end
372
373       subroutine xx2x(x,xx)
374       implicit real*8 (a-h,o-z)
375       include 'DIMENSIONS'
376       include 'COMMON.VAR'
377       include 'COMMON.CHAIN'
378       include 'COMMON.INTERACT'
379       double precision xx(maxvar),x(maxvar)
380
381       do i=1,nvar
382         x(i)=varall(i)
383       enddo
384
385       ig=0                                                                      
386       igall=0                                                                   
387       do i=4,nres                                                               
388         igall=igall+1                                                           
389         if (mask_phi(i).eq.1) then                                              
390           ig=ig+1                                                               
391           x(igall)=xx(ig)
392         endif                                                                   
393       enddo                                                                     
394                                                                                 
395       do i=3,nres                                                               
396         igall=igall+1                                                           
397         if (mask_theta(i).eq.1) then                                            
398           ig=ig+1                                                               
399           x(igall)=xx(ig)
400         endif                                                                   
401       enddo                                          
402
403       do ij=1,2                                                                 
404       do i=2,nres-1                                                             
405         if (itype(i).ne.10) then                                                
406           igall=igall+1                                                         
407           if (mask_side(i).eq.1) then                                           
408             ig=ig+1                                                             
409             x(igall)=xx(ig)
410           endif                                                                 
411         endif                                                                   
412       enddo                                                             
413       enddo                              
414
415       return
416       end
417   
418 c---------------------------------------------------------- 
419       subroutine minim_dc(etot,iretcode,nfun)
420       implicit real*8 (a-h,o-z)
421       include 'DIMENSIONS'
422       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
423 #ifdef MPI
424       include 'mpif.h'
425 #endif
426       include 'COMMON.SETUP'
427       include 'COMMON.IOUNITS'
428       include 'COMMON.VAR'
429       include 'COMMON.GEO'
430       include 'COMMON.MINIM'
431       include 'COMMON.CHAIN'
432       dimension iv(liv)                                               
433       double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
434 c      common /przechowalnia/ v
435
436       double precision energia(0:n_ene)
437       external func_dc,grad_dc,fdum
438       logical not_done,change,reduce 
439       double precision g(maxvar),f1
440
441       call deflt(2,iv,liv,lv,v)                                         
442 * 12 means fresh start, dont call deflt                                 
443       iv(1)=12                                                          
444 * max num of fun calls                                                  
445       if (maxfun.eq.0) maxfun=500
446       iv(17)=maxfun
447 * max num of iterations                                                 
448       if (maxmin.eq.0) maxmin=1000
449       iv(18)=maxmin
450 * controls output                                                       
451       iv(19)=2                                                          
452 * selects output unit                                                   
453       iv(21)=0
454       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout 
455 * 1 means to print out result                                           
456       iv(22)=print_min_res
457 * 1 means to print out summary stats                                    
458       iv(23)=print_min_stat
459 * 1 means to print initial x and d                                      
460       iv(24)=print_min_ini
461 * min val for v(radfac) default is 0.1                                  
462       v(24)=0.1D0                                                       
463 * max val for v(radfac) default is 4.0                                  
464       v(25)=2.0D0                                                       
465 c     v(25)=4.0D0                                                       
466 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
467 * the sumsl default is 0.1                                              
468       v(26)=0.1D0
469 * false conv if (act fnctn decrease) .lt. v(34)                         
470 * the sumsl default is 100*machep                                       
471       v(34)=v(34)/100.0D0                                               
472 * absolute convergence                                                  
473       if (tolf.eq.0.0D0) tolf=1.0D-4
474       v(31)=tolf
475 * relative convergence                                                  
476       if (rtolf.eq.0.0D0) rtolf=1.0D-4
477       v(32)=rtolf
478 * controls initial step size                                            
479        v(35)=1.0D-1                                                    
480 * large vals of d correspond to small components of step                
481       do i=1,6*nres
482         d(i)=1.0D-1
483       enddo
484
485       k=0
486       do i=1,nres-1
487         do j=1,3
488           k=k+1
489           x(k)=dc(j,i)
490         enddo
491       enddo
492       do i=2,nres-1
493         if (ialph(i,1).gt.0) then
494         do j=1,3
495           k=k+1
496           x(k)=dc(j,i+nres)
497         enddo
498         endif
499       enddo
500
501       call flush(iout)
502       call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)      
503
504       k=0
505       do i=1,nres-1
506         do j=1,3
507           k=k+1
508           dc(j,i)=x(k)
509         enddo
510       enddo
511       do i=2,nres-1
512         if (ialph(i,1).gt.0) then
513         do j=1,3
514           k=k+1
515           dc(j,i+nres)=x(k)
516         enddo
517         endif
518       enddo
519       call chainbuild_cart
520
521 cd      call zerograd
522 cd      nf=0
523 cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
524 cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
525 cd
526 cd      do i=1,k
527 cd       x(i)=x(i)+1.0D-5
528 cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
529 cd       x(i)=x(i)-1.0D-5
530 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
531 cd      enddo
532
533       etot=v(10)                                                      
534       iretcode=iv(1)
535       nfun=iv(6)
536       return  
537       end  
538
539       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)  
540       implicit real*8 (a-h,o-z)
541       include 'DIMENSIONS'
542 #ifdef MPI
543       include 'mpif.h'
544 #endif
545       include 'COMMON.SETUP'
546       include 'COMMON.DERIV'
547       include 'COMMON.IOUNITS'
548       include 'COMMON.GEO'
549       include 'COMMON.CHAIN'
550       include 'COMMON.VAR'
551       double precision energia(0:n_ene)
552       double precision ufparm
553       external ufparm                                                   
554       integer uiparm(1)                                                 
555       real*8 urparm(1)                                                    
556       dimension x(maxvar)
557       nfl=nf
558 cbad      icg=mod(nf,2)+1
559       icg=1
560
561       k=0
562       do i=1,nres-1
563         do j=1,3
564           k=k+1
565           dc(j,i)=x(k)
566         enddo
567       enddo
568       do i=2,nres-1
569         if (ialph(i,1).gt.0) then
570         do j=1,3
571           k=k+1
572           dc(j,i+nres)=x(k)
573         enddo
574         endif
575       enddo
576       call chainbuild_cart
577
578       call zerograd
579       call etotal(energia(0))
580       f=energia(0)
581
582 cd      print *,'func_dc ',nf,nfl,f
583
584       return                                                            
585       end                                                               
586
587       subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
588       implicit real*8 (a-h,o-z)
589       include 'DIMENSIONS'
590 #ifdef MPI
591       include 'mpif.h'
592 #endif
593       include 'COMMON.SETUP'
594       include 'COMMON.CHAIN'
595       include 'COMMON.DERIV'
596       include 'COMMON.VAR'
597       include 'COMMON.INTERACT'
598       include 'COMMON.FFIELD'
599       include 'COMMON.MD'
600       include 'COMMON.IOUNITS'
601       external ufparm
602       integer uiparm(1),k
603       double precision urparm(1)
604       dimension x(maxvar),g(maxvar)
605 c
606 c
607 c
608 cbad      icg=mod(nf,2)+1
609       icg=1
610 cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
611       if (nf-nfl+1) 20,30,40
612    20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
613 cd      print *,20
614       if (nf.eq.0) return
615       goto 40
616    30 continue
617 cd      print *,30
618       k=0
619       do i=1,nres-1
620         do j=1,3
621           k=k+1
622           dc(j,i)=x(k)
623         enddo
624       enddo
625       do i=2,nres-1
626         if (ialph(i,1).gt.0) then
627         do j=1,3
628           k=k+1
629           dc(j,i+nres)=x(k)
630         enddo
631         endif
632       enddo
633       call chainbuild_cart
634
635 C
636 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
637 C
638    40 call cartgrad
639 cd      print *,40
640       k=0
641       do i=1,nres-1
642         do j=1,3
643           k=k+1
644           g(k)=gcart(j,i)
645         enddo
646       enddo
647       do i=2,nres-1
648         if (ialph(i,1).gt.0) then
649         do j=1,3
650           k=k+1
651           g(k)=gxcart(j,i)
652         enddo
653         endif
654       enddo       
655
656       return
657       end