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
638 write (iout,*) "Variables set up nvarx",nvarx
639 write (iout,*) "Before energy minimization"
640 call etotal(energia(0))
641 call enerprint(energia(0))
646 c perform dynamic allocation of some global arrays
648 if (.not. allocated(scale)) allocate (scale(nvarx))
650 c set scaling parameter for function and derivative values;
651 c use square root of median eigenvalue of typical Hessian
663 c write (iout,*) "minim_dc Calling lbfgs"
664 call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
666 c write (iout,*) "minim_dc After lbfgs"
669 c write (iout,*) "checkgrad before SUMSL"
672 c call exec_checkgrad
675 call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
677 c write (iout,*) "checkgrad after SUMSL"
678 c call exec_checkgrad
689 if (ialph(i,1).gt.0) then
700 cd call func_dc(k,x,nf,f,idum,rdum,fdum)
701 cd call grad_dc(k,x,nf,g,idum,rdum,fdum)
705 cd call func_dc(k,x,nf,f1,idum,rdum,fdum)
707 cd print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
717 double precision function funcgrad_dc(x,g)
718 implicit real*8 (a-h,o-z)
723 include 'COMMON.SETUP'
724 include 'COMMON.CHAIN'
725 include 'COMMON.DERIV'
727 include 'COMMON.INTERACT'
728 include 'COMMON.FFIELD'
730 include 'COMMON.IOUNITS'
732 dimension x(maxvar),g(maxvar)
733 double precision energia(0:n_ene)
745 if (ialph(i,1).gt.0) then
754 call etotal(energia(0))
755 c write (iout,*) "energia",energia(0)
756 funcgrad_dc=energia(0)
758 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
769 if (ialph(i,1).gt.0) then
780 subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
781 implicit real*8 (a-h,o-z)
786 include 'COMMON.SETUP'
787 include 'COMMON.DERIV'
788 include 'COMMON.IOUNITS'
790 include 'COMMON.CHAIN'
792 double precision energia(0:n_ene)
793 double precision ufparm
810 if (ialph(i,1).gt.0) then
820 call etotal(energia(0))
823 cd print *,'func_dc ',nf,nfl,f
828 subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
829 implicit real*8 (a-h,o-z)
834 include 'COMMON.SETUP'
835 include 'COMMON.CHAIN'
836 include 'COMMON.DERIV'
838 include 'COMMON.INTERACT'
839 include 'COMMON.FFIELD'
841 include 'COMMON.IOUNITS'
844 double precision urparm(1)
845 dimension x(maxvar),g(maxvar)
851 cd print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
852 if (nf-nfl+1) 20,30,40
853 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
867 if (ialph(i,1).gt.0) then
877 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
889 if (ialph(i,1).gt.0) then