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