Merge branch 'bartek' of mmka.chem.univ.gda.pl:unres into bartek
[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       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           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           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           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           call fricmat_mult(z,d_a_tmp)
230         else if (iorder.eq.10) then
231           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)=0
439       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout 
440 * 1 means to print out result                                           
441       iv(22)=print_min_res
442 * 1 means to print out summary stats                                    
443       iv(23)=print_min_stat
444 * 1 means to print initial x and d                                      
445       iv(24)=print_min_ini
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