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:11)
201 c double precision Gcopy(maxres2,maxres2)
202 c common /przechowalnia/ Gcopy
204 integer i,j,iorder,ioverlap(maxres),ioverlap_last
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
304 else if (iorder.eq.11) then
305 call overlap_sc_list(ioverlap,ioverlap_last,.false.)
308 write (*,*) 'Processor',fg_rank,' CG group',kolor,
309 & ' absolute rank',myrank,' leves ERGASTULUM.'
310 write(*,*)'Processor',fg_rank,' wait times for respective orders',
311 & (' order[',i,']',time_order(i),i=0,10)
315 ************************************************************************
317 double precision function funcgrad(x,g)
320 include 'COMMON.CONTROL'
321 include 'COMMON.CHAIN'
322 include 'COMMON.DERIV'
324 include 'COMMON.INTERACT'
325 include 'COMMON.FFIELD'
327 include 'COMMON.QRESTR'
328 include 'COMMON.IOUNITS'
330 double precision energia(0:n_ene)
331 double precision x(nvar),g(nvar)
334 c write (iout,*) "in func x"
335 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
337 call var_to_geom(nvar,x)
339 call chainbuild_extconf
340 call etotal(energia(0))
343 call cart2intgrad(nvar,g)
345 C Add the components corresponding to local energy terms.
347 c Add the usampl contributions
350 gloc(i,icg)=gloc(i,icg)+dugamma(i)
353 gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
357 cd write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
358 g(i)=g(i)+gloc(i,icg)
363 subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
364 implicit real*8 (a-h,o-z)
366 include 'COMMON.DERIV'
367 include 'COMMON.IOUNITS'
370 double precision energia(0:n_ene)
372 double precision ufparm
378 c write (iout,*) "in func x"
379 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
383 cd print *,'func',nf,nfl,icg
384 call var_to_geom(n,x)
386 call chainbuild_extconf
387 cd write (iout,*) 'ETOTAL called from FUNC'
388 call etotal(energia(0))
392 c write (iout,*) "upon exit from func"
393 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
394 c write (iout,*) 'f=',f
399 ************************************************************************
400 subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
401 implicit real*8 (a-h,o-z)
403 include 'COMMON.DERIV'
404 include 'COMMON.IOUNITS'
407 double precision energia(0:n_ene)
409 double precision ufparm
415 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
419 call var_to_geom_restr(n,x)
421 call chainbuild_extconf
422 cd write (iout,*) 'ETOTAL called from FUNC'
423 call etotal(energia(0))
427 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
428 c write (iout,*) 'f=',etot
433 c-------------------------------------------------------
435 subroutine x2xx(x,xx,n)
436 implicit real*8 (a-h,o-z)
439 include 'COMMON.CHAIN'
440 include 'COMMON.INTERACT'
441 include 'COMMON.IOUNITS'
442 double precision xx(maxvar),x(maxvar)
444 c write (iout,*) "nvar",nvar
453 if (mask_phi(i).eq.1) then
461 if (mask_theta(i).eq.1) then
469 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
471 if (mask_side(i).eq.1) then
474 c write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
475 c write (iout,*) "x",x(igall)," xx",xx(ig)
486 subroutine xx2x(x,xx)
487 implicit real*8 (a-h,o-z)
490 include 'COMMON.CHAIN'
491 include 'COMMON.INTERACT'
492 include 'COMMON.IOUNITS'
493 double precision xx(maxvar),x(maxvar)
503 if (mask_phi(i).eq.1) then
511 if (mask_theta(i).eq.1) then
519 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
521 if (mask_side(i).eq.1) then
524 c write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
525 c write (iout,*) "x",x(igall)," xx",xx(ig)
534 c----------------------------------------------------------
535 subroutine minim_dc(etot,iretcode,nfun)
543 implicit real*8 (a-h,o-z)
546 parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
551 include 'COMMON.CONTROL'
552 include 'COMMON.SETUP'
553 include 'COMMON.IOUNITS'
556 include 'COMMON.MINIM'
557 include 'COMMON.CHAIN'
558 double precision minval,x(maxvar),d(maxvar),xx(maxvar)
560 double precision grdmin
561 double precision funcgrad_dc
562 external funcgrad_dc,optsave
565 double precision v(1:lv)
566 common /przechowalnia/ v
567 external func_dc,grad_dc,fdum
569 double precision g(maxvar),f1
571 double precision energia(0:n_ene)
574 coordtype='CARTESIAN'
577 jprint=print_min_stat
580 call deflt(2,iv,liv,lv,v)
581 * 12 means fresh start, dont call deflt
583 * max num of fun calls
584 if (maxfun.eq.0) maxfun=500
586 * max num of iterations
587 if (maxmin.eq.0) maxmin=1000
591 * selects output unit
593 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
594 * 1 means to print out result
596 * 1 means to print out summary stats
597 iv(23)=print_min_stat
598 * 1 means to print initial x and d
600 * min val for v(radfac) default is 0.1
602 * max val for v(radfac) default is 4.0
605 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
606 * the sumsl default is 0.1
608 * false conv if (act fnctn decrease) .lt. v(34)
609 * the sumsl default is 100*machep
611 * absolute convergence
612 if (tolf.eq.0.0D0) tolf=1.0D-4
614 * relative convergence
615 if (rtolf.eq.0.0D0) rtolf=1.0D-4
617 * controls initial step size
619 * large vals of d correspond to small components of step
632 if (ialph(i,1).gt.0) then
641 call etotal(energia(0))
642 call enerprint(energia(0))
647 c perform dynamic allocation of some global arrays
649 if (.not. allocated(scale)) allocate (scale(nvarx))
651 c set scaling parameter for function and derivative values;
652 c use square root of median eigenvalue of typical Hessian
664 c write (iout,*) "minim_dc Calling lbfgs"
665 call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
667 c write (iout,*) "minim_dc After lbfgs"
670 c write (iout,*) "checkgrad before SUMSL"
673 c call exec_checkgrad
676 call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
678 c write (iout,*) "checkgrad after SUMSL"
679 c call exec_checkgrad
690 if (ialph(i,1).gt.0) then
701 cd call func_dc(k,x,nf,f,idum,rdum,fdum)
702 cd call grad_dc(k,x,nf,g,idum,rdum,fdum)
706 cd call func_dc(k,x,nf,f1,idum,rdum,fdum)
708 cd print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
718 double precision function funcgrad_dc(x,g)
719 implicit real*8 (a-h,o-z)
724 include 'COMMON.SETUP'
725 include 'COMMON.CHAIN'
726 include 'COMMON.DERIV'
728 include 'COMMON.INTERACT'
729 include 'COMMON.FFIELD'
731 include 'COMMON.IOUNITS'
733 dimension x(maxvar),g(maxvar)
734 double precision energia(0:n_ene)
746 if (ialph(i,1).gt.0) then
755 call etotal(energia(0))
756 c write (iout,*) "energia",energia(0)
757 funcgrad_dc=energia(0)
759 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
770 if (ialph(i,1).gt.0) then
781 subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
782 implicit real*8 (a-h,o-z)
787 include 'COMMON.SETUP'
788 include 'COMMON.DERIV'
789 include 'COMMON.IOUNITS'
791 include 'COMMON.CHAIN'
793 double precision energia(0:n_ene)
794 double precision ufparm
811 if (ialph(i,1).gt.0) then
821 call etotal(energia(0))
824 cd print *,'func_dc ',nf,nfl,f
829 subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
830 implicit real*8 (a-h,o-z)
835 include 'COMMON.SETUP'
836 include 'COMMON.CHAIN'
837 include 'COMMON.DERIV'
839 include 'COMMON.INTERACT'
840 include 'COMMON.FFIELD'
842 include 'COMMON.IOUNITS'
845 double precision urparm(1)
846 dimension x(maxvar),g(maxvar)
852 cd print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
853 if (nf-nfl+1) 20,30,40
854 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
868 if (ialph(i,1).gt.0) then
878 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
890 if (ialph(i,1).gt.0) then