1 subroutine minimize(etot,x,iretcode,nfun)
5 parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
6 *********************************************************************
7 * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
8 * the calling subprogram. *
9 * when d(i)=1.0, then v(35) is the length of the initial step, *
10 * calculated in the usual pythagorean way. *
11 * absolute convergence occurs when the function is within v(31) of *
12 * zero. unless you know the minimum value in advance, abs convg *
13 * is probably not useful. *
14 * relative convergence is when the model predicts that the function *
15 * will decrease by less than v(32)*abs(fun). *
16 *********************************************************************
17 include 'COMMON.IOUNITS'
20 include 'COMMON.MINIM'
24 double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
25 double precision energia(0:n_ene)
29 external func,gradient,fdum
30 external func_restr,grad_restr
31 logical not_done,change,reduce
32 integer i,nvar_restr,nfun,iretcode
34 c common /przechowalnia/ v
42 call deflt(2,iv,liv,lv,v)
43 * 12 means fresh start, dont call deflt
45 * max num of fun calls
46 if (maxfun.eq.0) maxfun=500
48 * max num of iterations
49 if (maxmin.eq.0) maxmin=1000
55 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
56 * 1 means to print out result
58 * 1 means to print out summary stats
60 * 1 means to print initial x and d
62 * min val for v(radfac) default is 0.1
64 * max val for v(radfac) default is 4.0
67 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
68 * the sumsl default is 0.1
70 * false conv if (act fnctn decrease) .lt. v(34)
71 * the sumsl default is 100*machep
73 * absolute convergence
74 if (tolf.eq.0.0D0) tolf=1.0D-4
76 * relative convergence
77 if (rtolf.eq.0.0D0) rtolf=1.0D-4
79 * controls initial step size
81 * large vals of d correspond to small components of step
88 cd print *,'Calling SUMSL'
89 c call var_to_geom(nvar,x)
91 c call etotal(energia(0))
94 call x2xx(x,xx,nvar_restr)
95 call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
96 & iv,liv,lv,v,idum,rdum,fdum)
99 call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
103 cd print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
104 cd write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
107 call var_to_geom(nvar,x)
109 c write (iout,'(a)') 'Reduction worked, minimizing again...'
114 c call etotal(energia(0))
116 c call enerprint(energia(0))
119 c write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
126 c----------------------------------------------------------------------------
127 subroutine ergastulum
132 double precision time00
135 include 'COMMON.SETUP'
136 include 'COMMON.DERIV'
138 include 'COMMON.IOUNITS'
139 include 'COMMON.FFIELD'
140 include 'COMMON.INTERACT'
143 include 'COMMON.LAGRANGE.5diag'
145 include 'COMMON.LAGRANGE'
147 include 'COMMON.TIME1'
148 double precision z(maxres6),d_a_tmp(maxres6)
149 double precision edum(0:n_ene),time_order(0:10)
150 double precision Gcopy(maxres2,maxres2)
151 common /przechowalnia/ Gcopy
154 C Workers wait for variables and NF, and NFL from the boss
156 do while (iorder.ge.0)
157 c write (*,*) 'Processor',fg_rank,' CG group',kolor,
158 c & ' receives order from Master'
161 call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
162 time_Bcast=time_Bcast+MPI_Wtime()-time00
163 if (icall.gt.4 .and. iorder.ge.0)
164 & time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
168 c & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
169 if (iorder.eq.0) then
172 c write (2,*) "After etotal"
173 c write (2,*) "dimen",dimen," dimen3",dimen3
175 else if (iorder.eq.2) then
177 call etotal_short(edum)
178 c write (2,*) "After etotal_short"
179 c write (2,*) "dimen",dimen," dimen3",dimen3
181 else if (iorder.eq.3) then
183 call etotal_long(edum)
184 c write (2,*) "After etotal_long"
185 c write (2,*) "dimen",dimen," dimen3",dimen3
187 else if (iorder.eq.1) then
189 c write (2,*) "After sum_gradient"
190 c write (2,*) "dimen",dimen," dimen3",dimen3
193 else if (iorder.eq.4) then
194 call ginv_mult(z,d_a_tmp)
195 else if (iorder.eq.5) then
196 c Setup MD things for a slave
197 dimen=(nct-nnt+1)+nside
198 dimen1=(nct-nnt)+(nct-nnt+1)
200 c write (2,*) "dimen",dimen," dimen3",dimen3
202 call int_bounds(dimen,igmult_start,igmult_end)
203 igmult_start=igmult_start-1
204 call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
205 & ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
206 my_ng_count=igmult_end-igmult_start
207 call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
208 & MPI_INTEGER,FG_COMM,IERROR)
209 c write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
210 c write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
211 myginv_ng_count=maxres2*my_ng_count
212 c write (2,*) "igmult_start",igmult_start," igmult_end",
213 c & igmult_end," my_ng_count",my_ng_count
215 call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
216 & nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
217 call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
218 & nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
219 c write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
220 c write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
222 c call MPI_Barrier(FG_COMM,IERROR)
224 call MPI_Scatterv(ginv(1,1),nginv_counts(0),
225 & nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
226 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
228 time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
235 c write (2,*) "dimen",dimen," dimen3",dimen3
236 c write (2,*) "End MD setup"
238 c write (iout,*) "My chunk of ginv_block"
239 c call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
241 else if (iorder.eq.6) then
242 call int_from_cart1(.false.)
243 else if (iorder.eq.7) then
245 else if (iorder.eq.8) then
248 else if (iorder.eq.9) then
249 call fricmat_mult(z,d_a_tmp)
251 else if (iorder.eq.10) then
255 write (*,*) 'Processor',fg_rank,' CG group',kolor,
256 & ' absolute rank',myrank,' leves ERGASTULUM.'
257 write(*,*)'Processor',fg_rank,' wait times for respective orders',
258 & (' order[',i,']',time_order(i),i=0,10)
262 ************************************************************************
263 subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
264 implicit real*8 (a-h,o-z)
266 include 'COMMON.DERIV'
267 include 'COMMON.IOUNITS'
270 double precision energia(0:n_ene)
272 double precision ufparm
278 c write (iout,*) "in func x"
279 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
283 cd print *,'func',nf,nfl,icg
284 call var_to_geom(n,x)
286 call chainbuild_extconf
287 cd write (iout,*) 'ETOTAL called from FUNC'
288 call etotal(energia(0))
292 c write (iout,*) "upon exit from func"
293 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
294 c write (iout,*) 'f=',f
299 ************************************************************************
300 subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
301 implicit real*8 (a-h,o-z)
303 include 'COMMON.DERIV'
304 include 'COMMON.IOUNITS'
307 double precision energia(0:n_ene)
309 double precision ufparm
315 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
319 call var_to_geom_restr(n,x)
321 call chainbuild_extconf
322 cd write (iout,*) 'ETOTAL called from FUNC'
323 call etotal(energia(0))
327 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
328 c write (iout,*) 'f=',etot
333 c-------------------------------------------------------
334 subroutine x2xx(x,xx,n)
335 implicit real*8 (a-h,o-z)
338 include 'COMMON.CHAIN'
339 include 'COMMON.INTERACT'
340 double precision xx(maxvar),x(maxvar)
350 if (mask_phi(i).eq.1) then
358 if (mask_theta(i).eq.1) then
366 if (itype(i).ne.10) then
368 if (mask_side(i).eq.1) then
381 subroutine xx2x(x,xx)
382 implicit real*8 (a-h,o-z)
385 include 'COMMON.CHAIN'
386 include 'COMMON.INTERACT'
387 double precision xx(maxvar),x(maxvar)
397 if (mask_phi(i).eq.1) then
405 if (mask_theta(i).eq.1) then
413 if (itype(i).ne.10) then
415 if (mask_side(i).eq.1) then
426 c----------------------------------------------------------
427 subroutine minim_dc(etot,iretcode,nfun)
435 implicit real*8 (a-h,o-z)
438 parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
443 include 'COMMON.CONTROL'
444 include 'COMMON.SETUP'
445 include 'COMMON.IOUNITS'
448 include 'COMMON.MINIM'
449 include 'COMMON.CHAIN'
450 double precision minval,x(maxvar),d(maxvar),xx(maxvar)
452 double precision grdmin
453 double precision funcgrad_dc
457 double precision v(1:lv)
458 common /przechowalnia/ v
459 external func_dc,grad_dc,fdum
461 double precision g(maxvar),f1
465 coordtype='CARTESIAN'
468 jprint=print_min_stat
471 call deflt(2,iv,liv,lv,v)
472 * 12 means fresh start, dont call deflt
474 * max num of fun calls
475 if (maxfun.eq.0) maxfun=500
477 * max num of iterations
478 if (maxmin.eq.0) maxmin=1000
482 * selects output unit
484 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
485 * 1 means to print out result
487 * 1 means to print out summary stats
488 iv(23)=print_min_stat
489 * 1 means to print initial x and d
491 * min val for v(radfac) default is 0.1
493 * max val for v(radfac) default is 4.0
496 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
497 * the sumsl default is 0.1
499 * false conv if (act fnctn decrease) .lt. v(34)
500 * the sumsl default is 100*machep
502 * absolute convergence
503 if (tolf.eq.0.0D0) tolf=1.0D-4
505 * relative convergence
506 if (rtolf.eq.0.0D0) rtolf=1.0D-4
508 * controls initial step size
510 * large vals of d correspond to small components of step
523 if (ialph(i,1).gt.0) then
531 write (iout,*) "Variables set up nvarx",nvarx
536 c perform dynamic allocation of some global arrays
538 if (.not. allocated(scale)) allocate (scale(nvarx))
540 c set scaling parameter for function and derivative values;
541 c use square root of median eigenvalue of typical Hessian
553 write (iout,*) "Calling lbfgs"
554 call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
555 write (iout,*) "After lbfgs"
558 c write (iout,*) "checkgrad before SUMSL"
561 c call exec_checkgrad
564 call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
566 c write (iout,*) "checkgrad after SUMSL"
567 c call exec_checkgrad
578 if (ialph(i,1).gt.0) then
589 cd call func_dc(k,x,nf,f,idum,rdum,fdum)
590 cd call grad_dc(k,x,nf,g,idum,rdum,fdum)
594 cd call func_dc(k,x,nf,f1,idum,rdum,fdum)
596 cd print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
606 double precision function funcgrad_dc(x,g)
607 implicit real*8 (a-h,o-z)
612 include 'COMMON.SETUP'
613 include 'COMMON.CHAIN'
614 include 'COMMON.DERIV'
616 include 'COMMON.INTERACT'
617 include 'COMMON.FFIELD'
619 include 'COMMON.IOUNITS'
621 dimension x(maxvar),g(maxvar)
622 double precision energia(0:n_ene)
632 if (ialph(i,1).gt.0) then
641 call etotal(energia(0))
642 c write (iout,*) "energia",energia(0)
643 funcgrad_dc=energia(0)
645 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
657 if (ialph(i,1).gt.0) then
668 subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
669 implicit real*8 (a-h,o-z)
674 include 'COMMON.SETUP'
675 include 'COMMON.DERIV'
676 include 'COMMON.IOUNITS'
678 include 'COMMON.CHAIN'
680 double precision energia(0:n_ene)
681 double precision ufparm
698 if (ialph(i,1).gt.0) then
708 call etotal(energia(0))
711 cd print *,'func_dc ',nf,nfl,f
716 subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
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'
732 double precision urparm(1)
733 dimension x(maxvar),g(maxvar)
739 cd print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
740 if (nf-nfl+1) 20,30,40
741 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
755 if (ialph(i,1).gt.0) then
765 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
777 if (ialph(i,1).gt.0) then