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
571 coordtype='CARTESIAN'
574 jprint=print_min_stat
577 call deflt(2,iv,liv,lv,v)
578 * 12 means fresh start, dont call deflt
580 * max num of fun calls
581 if (maxfun.eq.0) maxfun=500
583 * max num of iterations
584 if (maxmin.eq.0) maxmin=1000
588 * selects output unit
590 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
591 * 1 means to print out result
593 * 1 means to print out summary stats
594 iv(23)=print_min_stat
595 * 1 means to print initial x and d
597 * min val for v(radfac) default is 0.1
599 * max val for v(radfac) default is 4.0
602 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
603 * the sumsl default is 0.1
605 * false conv if (act fnctn decrease) .lt. v(34)
606 * the sumsl default is 100*machep
608 * absolute convergence
609 if (tolf.eq.0.0D0) tolf=1.0D-4
611 * relative convergence
612 if (rtolf.eq.0.0D0) rtolf=1.0D-4
614 * controls initial step size
616 * large vals of d correspond to small components of step
629 if (ialph(i,1).gt.0) then
637 write (iout,*) "Variables set up nvarx",nvarx
642 c perform dynamic allocation of some global arrays
644 if (.not. allocated(scale)) allocate (scale(nvarx))
646 c set scaling parameter for function and derivative values;
647 c use square root of median eigenvalue of typical Hessian
659 c write (iout,*) "minim_dc Calling lbfgs"
660 call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
662 c write (iout,*) "minim_dc After lbfgs"
665 c write (iout,*) "checkgrad before SUMSL"
668 c call exec_checkgrad
671 call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
673 c write (iout,*) "checkgrad after SUMSL"
674 c call exec_checkgrad
685 if (ialph(i,1).gt.0) then
696 cd call func_dc(k,x,nf,f,idum,rdum,fdum)
697 cd call grad_dc(k,x,nf,g,idum,rdum,fdum)
701 cd call func_dc(k,x,nf,f1,idum,rdum,fdum)
703 cd print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
713 double precision function funcgrad_dc(x,g)
714 implicit real*8 (a-h,o-z)
719 include 'COMMON.SETUP'
720 include 'COMMON.CHAIN'
721 include 'COMMON.DERIV'
723 include 'COMMON.INTERACT'
724 include 'COMMON.FFIELD'
726 include 'COMMON.IOUNITS'
728 dimension x(maxvar),g(maxvar)
729 double precision energia(0:n_ene)
741 if (ialph(i,1).gt.0) then
750 call etotal(energia(0))
751 c write (iout,*) "energia",energia(0)
752 funcgrad_dc=energia(0)
754 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
765 if (ialph(i,1).gt.0) then
776 subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
777 implicit real*8 (a-h,o-z)
782 include 'COMMON.SETUP'
783 include 'COMMON.DERIV'
784 include 'COMMON.IOUNITS'
786 include 'COMMON.CHAIN'
788 double precision energia(0:n_ene)
789 double precision ufparm
806 if (ialph(i,1).gt.0) then
816 call etotal(energia(0))
819 cd print *,'func_dc ',nf,nfl,f
824 subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
825 implicit real*8 (a-h,o-z)
830 include 'COMMON.SETUP'
831 include 'COMMON.CHAIN'
832 include 'COMMON.DERIV'
834 include 'COMMON.INTERACT'
835 include 'COMMON.FFIELD'
837 include 'COMMON.IOUNITS'
840 double precision urparm(1)
841 dimension x(maxvar),g(maxvar)
847 cd print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
848 if (nf-nfl+1) 20,30,40
849 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
863 if (ialph(i,1).gt.0) then
873 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
885 if (ialph(i,1).gt.0) then