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