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