1 subroutine minimize(etot,x,iretcode,nfun)
13 parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
15 *********************************************************************
16 * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
17 * the calling subprogram. *
18 * when d(i)=1.0, then v(35) is the length of the initial step, *
19 * calculated in the usual pythagorean way. *
20 * absolute convergence occurs when the function is within v(31) of *
21 * zero. unless you know the minimum value in advance, abs convg *
22 * is probably not useful. *
23 * relative convergence is when the model predicts that the function *
24 * will decrease by less than v(32)*abs(fun). *
25 *********************************************************************
26 include 'COMMON.IOUNITS'
29 include 'COMMON.MINIM'
33 double precision grdmin
38 double precision v(1:lv)
39 common /przechowalnia/ v
43 external func,gradient,fdum
44 external func_restr,grad_restr
45 logical not_done,change,reduce
47 double precision x(maxvar),d(maxvar),xx(maxvar)
48 double precision energia(0:n_ene)
49 integer i,nvar_restr,nfun,iretcode
51 c common /przechowalnia/ v
60 if (.not. allocated(scale)) allocate (scale(nvar))
62 c set scaling parameter for function and derivative values;
63 c use square root of median eigenvalue of typical Hessian
75 c write (iout,*) "Calling lbfgs"
76 write (iout,*) 'Calling LBFGS, minimization in angles'
77 call var_to_geom(nvar,x)
78 call chainbuild_extconf
79 call etotal(energia(0))
80 call enerprint(energia(0))
81 call lbfgs (nvar,x,etot,grdmin,funcgrad,optsave)
83 write (iout,*) "Minimized energy",etot
91 call deflt(2,iv,liv,lv,v)
92 * 12 means fresh start, dont call deflt
94 * max num of fun calls
95 if (maxfun.eq.0) maxfun=500
97 * max num of iterations
98 if (maxmin.eq.0) maxmin=1000
102 * selects output unit
104 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
105 * 1 means to print out result
107 * 1 means to print out summary stats
108 iv(23)=print_min_stat
109 * 1 means to print initial x and d
111 * min val for v(radfac) default is 0.1
113 * max val for v(radfac) default is 4.0
116 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
117 * the sumsl default is 0.1
119 * false conv if (act fnctn decrease) .lt. v(34)
120 * the sumsl default is 100*machep
122 * absolute convergence
123 if (tolf.eq.0.0D0) tolf=1.0D-4
125 * relative convergence
126 if (rtolf.eq.0.0D0) rtolf=1.0D-4
128 * controls initial step size
130 * large vals of d correspond to small components of step
137 write (iout,*) 'Calling SUMSL'
138 call var_to_geom(nvar,x)
139 call chainbuild_extconf
141 call etotal(energia(0))
142 call enerprint(energia(0))
145 call x2xx(x,xx,nvar_restr)
146 call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
147 & iv,liv,lv,v,idum,rdum,fdum)
150 call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
154 cd print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
155 cd write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
158 call var_to_geom(nvar,x)
160 c write (iout,'(a)') 'Reduction worked, minimizing again...'
164 call chainbuild_extconf
165 c call etotal(energia(0))
167 c call enerprint(energia(0))
170 c write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
177 c----------------------------------------------------------------------------
178 subroutine ergastulum
183 double precision time00
186 include 'COMMON.SETUP'
187 include 'COMMON.DERIV'
189 include 'COMMON.IOUNITS'
190 include 'COMMON.FFIELD'
191 include 'COMMON.INTERACT'
194 include 'COMMON.LAGRANGE.5diag'
196 include 'COMMON.LAGRANGE'
198 include 'COMMON.TIME1'
199 double precision z(maxres6),d_a_tmp(maxres6)
200 double precision edum(0:n_ene),time_order(0:10)
201 double precision Gcopy(maxres2,maxres2)
202 common /przechowalnia/ Gcopy
205 C Workers wait for variables and NF, and NFL from the boss
207 do while (iorder.ge.0)
208 c write (*,*) 'Processor',fg_rank,' CG group',kolor,
209 c & ' receives order from Master'
212 call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
213 time_Bcast=time_Bcast+MPI_Wtime()-time00
214 if (icall.gt.4 .and. iorder.ge.0)
215 & time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
219 c & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
220 if (iorder.eq.0) then
223 c write (2,*) "After etotal"
224 c write (2,*) "dimen",dimen," dimen3",dimen3
226 else if (iorder.eq.2) then
228 call etotal_short(edum)
229 c write (2,*) "After etotal_short"
230 c write (2,*) "dimen",dimen," dimen3",dimen3
232 else if (iorder.eq.3) then
234 call etotal_long(edum)
235 c write (2,*) "After etotal_long"
236 c write (2,*) "dimen",dimen," dimen3",dimen3
238 else if (iorder.eq.1) then
240 c write (2,*) "After sum_gradient"
241 c write (2,*) "dimen",dimen," dimen3",dimen3
244 else if (iorder.eq.4) then
245 call ginv_mult(z,d_a_tmp)
246 else if (iorder.eq.5) then
247 c Setup MD things for a slave
248 dimen=(nct-nnt+1)+nside
249 dimen1=(nct-nnt)+(nct-nnt+1)
251 c write (2,*) "dimen",dimen," dimen3",dimen3
253 call int_bounds(dimen,igmult_start,igmult_end)
254 igmult_start=igmult_start-1
255 call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
256 & ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
257 my_ng_count=igmult_end-igmult_start
258 call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
259 & MPI_INTEGER,FG_COMM,IERROR)
260 c write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
261 c write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
262 myginv_ng_count=maxres2*my_ng_count
263 c write (2,*) "igmult_start",igmult_start," igmult_end",
264 c & igmult_end," my_ng_count",my_ng_count
266 call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
267 & nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
268 call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
269 & nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
270 c write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
271 c write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
273 c call MPI_Barrier(FG_COMM,IERROR)
275 call MPI_Scatterv(ginv(1,1),nginv_counts(0),
276 & nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
277 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
279 time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
286 c write (2,*) "dimen",dimen," dimen3",dimen3
287 c write (2,*) "End MD setup"
289 c write (iout,*) "My chunk of ginv_block"
290 c call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
292 else if (iorder.eq.6) then
293 call int_from_cart1(.false.)
294 else if (iorder.eq.7) then
296 else if (iorder.eq.8) then
299 else if (iorder.eq.9) then
300 call fricmat_mult(z,d_a_tmp)
302 else if (iorder.eq.10) then
306 write (*,*) 'Processor',fg_rank,' CG group',kolor,
307 & ' absolute rank',myrank,' leves ERGASTULUM.'
308 write(*,*)'Processor',fg_rank,' wait times for respective orders',
309 & (' order[',i,']',time_order(i),i=0,10)
313 ************************************************************************
315 double precision function funcgrad(x,g)
318 include 'COMMON.CONTROL'
319 include 'COMMON.CHAIN'
320 include 'COMMON.DERIV'
322 include 'COMMON.INTERACT'
323 include 'COMMON.FFIELD'
325 include 'COMMON.QRESTR'
326 include 'COMMON.IOUNITS'
328 double precision energia(0:n_ene)
329 double precision x(nvar),g(nvar)
332 c write (iout,*) "in func x"
333 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
335 call var_to_geom(nvar,x)
337 call chainbuild_extconf
338 call etotal(energia(0))
341 call cart2intgrad(nvar,g)
343 C Add the components corresponding to local energy terms.
345 c Add the usampl contributions
348 gloc(i,icg)=gloc(i,icg)+dugamma(i)
351 gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
355 cd write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
356 g(i)=g(i)+gloc(i,icg)
361 subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
362 implicit real*8 (a-h,o-z)
364 include 'COMMON.DERIV'
365 include 'COMMON.IOUNITS'
368 double precision energia(0:n_ene)
370 double precision ufparm
376 c write (iout,*) "in func x"
377 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
381 cd print *,'func',nf,nfl,icg
382 call var_to_geom(n,x)
384 call chainbuild_extconf
385 cd write (iout,*) 'ETOTAL called from FUNC'
386 call etotal(energia(0))
390 c write (iout,*) "upon exit from func"
391 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
392 c write (iout,*) 'f=',f
397 ************************************************************************
398 subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
399 implicit real*8 (a-h,o-z)
401 include 'COMMON.DERIV'
402 include 'COMMON.IOUNITS'
405 double precision energia(0:n_ene)
407 double precision ufparm
413 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
417 call var_to_geom_restr(n,x)
419 call chainbuild_extconf
420 cd write (iout,*) 'ETOTAL called from FUNC'
421 call etotal(energia(0))
425 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
426 c write (iout,*) 'f=',etot
431 c-------------------------------------------------------
433 subroutine x2xx(x,xx,n)
434 implicit real*8 (a-h,o-z)
437 include 'COMMON.CHAIN'
438 include 'COMMON.INTERACT'
439 include 'COMMON.IOUNITS'
440 double precision xx(maxvar),x(maxvar)
442 c write (iout,*) "nvar",nvar
451 if (mask_phi(i).eq.1) then
459 if (mask_theta(i).eq.1) then
467 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
469 if (mask_side(i).eq.1) then
472 c write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
473 c write (iout,*) "x",x(igall)," xx",xx(ig)
484 subroutine xx2x(x,xx)
485 implicit real*8 (a-h,o-z)
488 include 'COMMON.CHAIN'
489 include 'COMMON.INTERACT'
490 include 'COMMON.IOUNITS'
491 double precision xx(maxvar),x(maxvar)
501 if (mask_phi(i).eq.1) then
509 if (mask_theta(i).eq.1) then
517 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
519 if (mask_side(i).eq.1) then
522 c write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
523 c write (iout,*) "x",x(igall)," xx",xx(ig)
532 c----------------------------------------------------------
533 subroutine minim_dc(etot,iretcode,nfun)
541 implicit real*8 (a-h,o-z)
544 parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
549 include 'COMMON.CONTROL'
550 include 'COMMON.SETUP'
551 include 'COMMON.IOUNITS'
554 include 'COMMON.MINIM'
555 include 'COMMON.CHAIN'
556 double precision minval,x(maxvar),d(maxvar),xx(maxvar)
558 double precision grdmin
559 double precision funcgrad_dc
560 external funcgrad_dc,optsave
563 double precision v(1:lv)
564 common /przechowalnia/ v
565 external func_dc,grad_dc,fdum
567 double precision g(maxvar),f1
569 double precision energia(0:n_ene)
572 coordtype='CARTESIAN'
575 jprint=print_min_stat
578 call deflt(2,iv,liv,lv,v)
579 * 12 means fresh start, dont call deflt
581 * max num of fun calls
582 if (maxfun.eq.0) maxfun=500
584 * max num of iterations
585 if (maxmin.eq.0) maxmin=1000
589 * selects output unit
591 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
592 * 1 means to print out result
594 * 1 means to print out summary stats
595 iv(23)=print_min_stat
596 * 1 means to print initial x and d
598 * min val for v(radfac) default is 0.1
600 * max val for v(radfac) default is 4.0
603 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
604 * the sumsl default is 0.1
606 * false conv if (act fnctn decrease) .lt. v(34)
607 * the sumsl default is 100*machep
609 * absolute convergence
610 if (tolf.eq.0.0D0) tolf=1.0D-4
612 * relative convergence
613 if (rtolf.eq.0.0D0) rtolf=1.0D-4
615 * controls initial step size
617 * large vals of d correspond to small components of step
630 if (ialph(i,1).gt.0) then
639 call etotal(energia(0))
640 call enerprint(energia(0))
645 c perform dynamic allocation of some global arrays
647 if (.not. allocated(scale)) allocate (scale(nvarx))
649 c set scaling parameter for function and derivative values;
650 c use square root of median eigenvalue of typical Hessian
662 c write (iout,*) "minim_dc Calling lbfgs"
663 call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
665 c write (iout,*) "minim_dc After lbfgs"
668 c write (iout,*) "checkgrad before SUMSL"
671 c call exec_checkgrad
674 call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
676 c write (iout,*) "checkgrad after SUMSL"
677 c call exec_checkgrad
688 if (ialph(i,1).gt.0) then
699 cd call func_dc(k,x,nf,f,idum,rdum,fdum)
700 cd call grad_dc(k,x,nf,g,idum,rdum,fdum)
704 cd call func_dc(k,x,nf,f1,idum,rdum,fdum)
706 cd print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
716 double precision function funcgrad_dc(x,g)
717 implicit real*8 (a-h,o-z)
722 include 'COMMON.SETUP'
723 include 'COMMON.CHAIN'
724 include 'COMMON.DERIV'
726 include 'COMMON.INTERACT'
727 include 'COMMON.FFIELD'
729 include 'COMMON.IOUNITS'
731 dimension x(maxvar),g(maxvar)
732 double precision energia(0:n_ene)
744 if (ialph(i,1).gt.0) then
753 call etotal(energia(0))
754 c write (iout,*) "energia",energia(0)
755 funcgrad_dc=energia(0)
757 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
768 if (ialph(i,1).gt.0) then
779 subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
780 implicit real*8 (a-h,o-z)
785 include 'COMMON.SETUP'
786 include 'COMMON.DERIV'
787 include 'COMMON.IOUNITS'
789 include 'COMMON.CHAIN'
791 double precision energia(0:n_ene)
792 double precision ufparm
809 if (ialph(i,1).gt.0) then
819 call etotal(energia(0))
822 cd print *,'func_dc ',nf,nfl,f
827 subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
828 implicit real*8 (a-h,o-z)
833 include 'COMMON.SETUP'
834 include 'COMMON.CHAIN'
835 include 'COMMON.DERIV'
837 include 'COMMON.INTERACT'
838 include 'COMMON.FFIELD'
840 include 'COMMON.IOUNITS'
843 double precision urparm(1)
844 dimension x(maxvar),g(maxvar)
850 cd print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
851 if (nf-nfl+1) 20,30,40
852 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
866 if (ialph(i,1).gt.0) then
876 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
888 if (ialph(i,1).gt.0) then