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