1 subroutine maxlikopt(nvarr,x,xmin,f0)
4 include "DIMENSIONS.ZSCOPT"
10 include "COMMON.WEIGHTS"
11 include "COMMON.VMCPAR"
12 include "COMMON.OPTIM"
13 include "COMMON.CLASSES"
14 include "COMMON.ENERGIES"
15 include "COMMON.IOUNITS"
16 integer i,j,k,iprot,nf,ii,inn,n,nn,indn(max_ene)
17 logical solved,all_satisfied,not_done
18 double precision ran_number,f,f0,x(max_paropt),xmin(max_paropt),
21 integer iran_num,nvarr
22 double precision tcpu,t0,t1,t0w,t1w
29 if (me.eq.master) then
30 if (nparmset.gt.0) then
31 write (liczba,'(bz,i3.3)') myparm
32 restartname=prefix(:lenpre)//"_par"//liczba//'.restart'
34 restartname=prefix(:lenpre)//'.restart'
38 restartname=prefix(:lenpre)//'.restart'
40 print *,"Start solving inequalities"
43 c Only evaluate initial energies
48 call minimize_setup(nvarr,x)
51 write (iout,*) "CPU time for setup:",t1-t0," sec."
53 write (iout,*) "Wall clock time for setup:",t1w-t0w," sec."
60 call print_minimize_results(nvarr)
68 call minimize_setup(nvarr,x)
71 write (iout,*) "CPU time for setup:",t1-t0," sec."
73 write (iout,*) "Wall clock time for setup:",t1w-t0w," sec."
77 print *,"Calling minimize"
81 call minimize(nvarr,xmin,f0)
83 write (iout,*) "CPU time for minimization:",t1-t0," sec."
86 write (iout,*) "Wall clock time for minimization:",t1w-t0w
93 c write (iout,*) "Start MCM procedure"
96 call print_minimize_results(nvarr)
103 print *,"Calling scan"
104 call minimize_setup(nvarr,x)
108 call scan(nvarr,xmin,f0,.true.)
109 call print_minimize_results(nvarr)
112 60 format(bz,100f10.5)
115 c-----------------------------------------------------------------------------
116 subroutine minimize_setup(nvar,x)
119 include "DIMENSIONS.ZSCOPT"
122 integer IERROR,TAG,STATUS(MPI_STATUS_SIZE)
125 include "COMMON.IOUNITS"
126 include "COMMON.CLASSES"
127 include "COMMON.NAMES"
128 include "COMMON.OPTIM"
129 include "COMMON.WEIGHTS"
130 include "COMMON.WEIGHTDER"
131 include "COMMON.VMCPAR"
132 include "COMMON.ENERGIES"
133 include "COMMON.TIME1"
134 double precision viol
135 integer i,j,k,ii,ib,ibatch,iprot,nvar,nf
136 double precision x(max_paropt)
140 write (iout,*) "Entered MINIMIZE_SETUP"
143 call func1(nvar,3,x_orig)
148 c write (iout,*) "Initial X",(x(i),i=1,nvar)
149 c write (iout,*) "X_orig",(x_orig(i),i=1,nvar)
150 c write (iout,*) "ww_oorig",(ww_oorig(i),i=1,n_ene)
154 write (iout,*) "x",(x(i),i=1,nvar)
157 write(iout,'(10(1x,a6,1x))')(wname(print_order(i))(:6),i=1,n_ene)
158 write(iout,40)(ww(print_order(i)),i=1,n_ene)
159 write(iout,*)'-----------------------------------'
163 40 format(10(f7.4,1x))
165 c------------------------------------------------------------------------------
166 subroutine print_minimize_results(nvarr)
169 include "DIMENSIONS.ZSCOPT"
172 integer ierror,status(MPI_STATUS_SIZE),tag
175 include "COMMON.CONTROL"
176 include "COMMON.WEIGHTS"
177 include "COMMON.WEIGHTDER"
178 include "COMMON.PROTNAME"
179 include "COMMON.ENERGIES"
180 include "COMMON.CLASSES"
181 include "COMMON.IOUNITS"
182 include "COMMON.VMCPAR"
183 include "COMMON.OPTIM"
184 include "COMMON.TIME1"
185 integer nf,nvarr,i,ib,ic,ii,num_nat_tot,iprot,igather
186 double precision x(max_ene)
187 double precision viol
191 c write (iout,*) "print_minimize_results"
193 write(iout,*) '----------------------'
195 call x2w(nvarr,xm(1))
196 c write (iout,*) "Calling func1"
198 call func1(nvarr,2,xm(1))
199 c write (iout,*) "Calling write_params"
201 call write_params(nvarr,nf,xm(1))
202 c write(iout,*) "After write_params"
206 if (nparmset.eq.1) then
207 call write_prot(ii,'NewE_all')
209 write (liczba,'(bz,i3.3)') myparm
210 call write_prot(ii,'NewE_all'//liczba)
213 c call make_distrib(iclass(ic,ii),ii)
221 c------------------------------------------------------------------------------
222 subroutine write_params(nvar,nf,x)
225 include "DIMENSIONS.ZSCOPT"
226 include "COMMON.CONTROL"
227 include "COMMON.WEIGHTS"
228 include "COMMON.NAMES"
229 include "COMMON.INTERACT"
230 include "COMMON.TORSION"
231 include "COMMON.SCCOR"
232 include "COMMON.CLASSES"
233 include "COMMON.ENERGIES"
234 include "COMMON.IOUNITS"
235 include "COMMON.FFIELD"
236 include "COMMON.OPTIM"
238 double precision x(max_paropt)
241 integer i,ii,j,l,ll,jj,pj,jstart,jend,it1,it2,k,iblock
242 integer itor2typ(-2:2) /-20,9,10,9,20/
243 character*1 aster(0:1) /" ","*"/,ast11,ast12,ast21,ast22
245 character*1 toronelet(-2:2) /"p","a","G","A","P"/
249 write (iout,*) "Weights (asterisk: optimized)"
250 write(iout,'(1x,15(a8,2x))') (wname(print_order(i))(:8),i=1,n_ene)
251 write(iout,50)(ww(print_order(i)),
252 & aster(imask(print_order(i))),i=1,n_ene)
254 C Write weights in UNRES format for the input file
255 if (nparmset.eq.1) then
256 open(88,file="weights_opt."//prefix(:ilen(prefix)),
259 write (liczba,'(bz,i3.3)') myparm
260 open(88,file="weights_opt."//prefix(:ilen(prefix))//"par_"
270 write(88,'(bz,2a,f7.5," ",$)')
271 & wname(pj)(:ilen(wname(pj))),'=',ww(pj)
272 jj=jj+ilen(wname(pj))+9
274 write (88,'(a,$)') (" ",j=1,79-jj)
277 write (88,'(bz,2(2a,f7.5," "))')
278 & "CUTOFF","=",7.0d0,"WCORR4","=",0.0d0
280 if (mod_side .or. mod_side_other) then
281 if (nparmset.eq.1) then
283 & file="sc_"//POT(:ilen(POT))//"_opt."//prefix(:ilen(prefix)),
287 & file="sc_"//POT(:ilen(POT))//"_opt."//prefix(:ilen(prefix))//
292 write(88,'(2i5)') ipot,expon
293 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
294 & ', exponents are ',expon,2*expon
295 goto (10,20,30,30,40) ipot
296 C----------------------- LJ potential ---------------------------------
297 10 write (88,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
299 write (iout,'(/a/)') 'Parameters of the LJ potential:'
300 write (iout,'(a/)') 'The epsilon array:'
301 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
302 write (iout,'(/a)') 'One-body parameters:'
303 write (iout,'(a,4x,a)') 'residue','sigma'
304 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
307 C----------------------- LJK potential --------------------------------
308 20 write (88,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
309 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
311 write (iout,'(/a/)') 'Parameters of the LJK potential:'
312 write (iout,'(a/)') 'The epsilon array:'
313 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
314 write (iout,'(/a)') 'One-body parameters:'
315 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
316 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
320 C---------------------- GB or BP potential -----------------------------
321 30 write (88,'(4f20.15)')((eps(i,j),j=i,ntyp),i=1,ntyp),
322 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
324 C For the GB potential convert sigma'**2 into chi'
327 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
332 write (iout,'(/a/)') 'Parameters of the BP potential:'
334 write (iout,'(/a/)') 'Parameters of the GB potential:'
336 write (iout,'(a/)') 'The epsilon array:'
337 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
338 write (iout,'(/a)') 'One-body parameters (asterisk: optimized):'
339 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
342 write (iout,'(a3,6x,f10.5,1x,a1,3(f8.5,1x,a1))') restyp(i),
343 & sigma0(i),aster(mask_sigma(1,i)),sigii(i),aster(mask_sigma(2,i)),
344 & chip(i),aster(mask_sigma(3,i)),alp(i),aster(mask_sigma(4,i))
348 C--------------------- GBV potential -----------------------------------
349 40 write (88,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
350 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
351 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
353 write (iout,'(/a/)') 'Parameters of the GBV potential:'
354 write (iout,'(a/)') 'The epsilon array:'
355 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
356 write (iout,'(/a)') 'One-body parameters:'
357 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
358 & 's||/s_|_^2',' chip ',' alph '
359 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
360 & sigii(i),chip(i),alp(i),i=1,ntyp)
363 C-----------------------------------------------------------------------
367 if (mod_tor .or. tor_mode.eq.2) then
371 write (iout,'(a)') "Optimized weights of the torsional terms."
373 do i=-ntortyp+1,ntortyp-1
374 do j=-ntortyp+1,ntortyp-1
375 if (mask_tor(0,j,i,iblock).gt.0) then
376 write (iout,'(i4,2a4,f10.5)') 0,restyp(itor2typ(j)),
377 & restyp(itor2typ(i)),weitor(0,j,i,iblock)
382 if (nparmset.eq.1) then
384 open(88,file="tor_opt.parm."//prefix(:ilen(prefix)),
386 write (88,'(i1,7h *** ,2a,7h *** )') ntortyp,
387 & "Optimized torsional parameters ",prefix(:ilen(prefix))
388 write (88,'(40i2)') (itortyp(i),i=1,ntyp)
390 if (tor_mode.eq.0) then
392 write (iout,'(/a/)') 'Torsional constants:'
396 write (iout,*) 'ityp',i,' jtyp',j
397 write (iout,*) 'Fourier constants'
398 do k=1,nterm(i,j,iblock)
399 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
402 write (iout,*) 'Lorenz constants'
403 do k=1,nlor(i,j,iblock)
404 write (iout,'(3(1pe15.5))')
405 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
413 do j=-ntortyp+1,ntortyp-1
414 write (88,'(2i2,16h ************ ,a1,1h-,a1)')
415 & nterm(i,j,iblock),nlor(i,j,iblock),onelet(itor2typ(i)),
416 & onelet(itor2typ(j))
417 do k=1,nterm(i,j,iblock)
418 write (88,'(i6,13x,2(1pe18.5))') k,
419 & v1(k,i,j,iblock)*weitor(0,i,j,iblock),v2(k,i,j,
420 & iblock)*weitor(0,i,j,iblock)
422 do k=1,nlor(i,j,iblock)
423 write (88,'(i6,13x,3(1pe18.5))')
424 & k,vlor1(k,i,j)*weitor(k,i,j,iblock),
425 & vlor2(k,i,j)*weitor(k,i,j,iblock),
426 & vlor3(k,i,j)*weitor(k,i,j,iblock)
433 c Print valence-torsional parameters
435 & "Parameters of the valence-torsional potentials"
436 do i=-ntortyp+1,ntortyp-1
437 do j=-ntortyp+1,ntortyp-1
438 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
439 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
440 do k=1,nterm_kcc(j,i)
441 do l=1,nterm_kcc_Tb(j,i)
442 do ll=1,nterm_kcc_Tb(j,i)
443 write (iout,'(3i5,2f15.4)')
444 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
451 do i=-ntortyp+1,ntortyp-1
452 do j=-ntortyp+1,ntortyp-1
453 write (88,'(16h ************ ,a1,1h-,a1)')
454 & onelet(itor2typ(i)),onelet(itor2typ(j))
455 write (88,'(2i5)') nterm_kcc(j,i),nterm_kcc_Tb(j,i)
456 do k=1,nterm_kcc(j,i)
457 do l=1,nterm_kcc_Tb(j,i)
458 do ll=1,nterm_kcc_Tb(j,i)
459 write (88,'(3i5,2(1pe15.5))') k,l-1,ll-1,
460 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
474 write (iout,'(a)') "Optimized weights of the sccor terms."
478 if (mask_tor(l,j,i,iblock).gt.0) then
479 write (iout,'(i4,2a4,f10.5)') l,restyp(isccortyp(j)),
480 & restyp(isccortyp(i)),weitor(l,j,i,1)
486 if (nparmset.eq.1) then
488 open(88,file="sccor_opt.parm."//prefix(:ilen(prefix)),
490 write (88,'(i2,7h *** ,2a,7h *** )') nsccortyp,
491 & "Optimized local correlation parameters ",prefix(:ilen(prefix))
492 write (88,'(20i4)') (isccortyp(i),i=1,ntyp)
496 write (88,'(2i4)') nterm_sccor(i,j),nlor_sccor(i,j)
497 write (88,'(12(1h*),1x,a1,1h-,a1)') onelet(i),onelet(j)
498 do k=1,nterm_sccor(i,j)
499 write (88,'(i6,13x,2(1pe18.5))') k,
500 & v1sccor(k,l,i,j)*weitor(l,i,j,1),
501 & v2sccor(k,l,i,j)*weitor(l,i,j,1)
511 C-----------------------------------------------------------------------
512 if (mod_fourier(nloctyp).gt.0) then
514 & "Fourier coefficients of cumulants (asterisk: optimized)"
517 write (iout,*) "Type: ",onelet(iloctyp(i))
518 write (iout,*) "Coefficients of the expansion of B1"
520 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
522 write (iout,*) "Coefficients of the expansion of B2"
524 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
526 write (iout,*) "Coefficients of the expansion of C"
527 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
528 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
529 write (iout,*) "Coefficients of the expansion of D"
530 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
531 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
532 write (iout,*) "Coefficients of the expansion of E"
533 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
536 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
540 write (iout,*) 'Type ',restyp(iloctyp(i))
541 write (iout,'(a,i2,a,f10.5,1x,a)')
542 & ('b(',j,')=',b(j,i),aster(mask_fourier(j,i)),j=1,13)
545 c Write a file with parameters for UNRES
546 if (nparmset.eq.1) then
547 open(88,file="fourier_opt.parm."//prefix(:ilen(prefix)),
550 open(88,file="fourier_opt.parm."//prefix(:ilen(prefix))//
554 write (88,'(i5,5x,a)') nloctyp,
555 & "# Number of local interaction types"
557 write (88,'(40i2)') (itype2loc(i),i=1,ntyp)
558 write (88,'(20i4)') (iloctyp(i),i=0,nloctyp)
560 write (88,'(a)') restyp(iloctyp(i))
563 write (88,'(f20.10,4h b1,i1,1h(,i1,1h))')
569 write (88,'(f20.10,4h b2,i1,1h(,i1,1h))')
575 write (88,'(f20.10,4h c1,i1,1h(,i1,1h))')
576 & 2*ccnew(ii,j,i),j,ii
581 write (88,'(f20.10,4h d1,i1,1h(,i1,1h))')
582 & 2*ddnew(ii,j,i),j,ii
588 write (88,'(f20.10,3h e,2i1,1h(,i1,1h))')
589 & eenew(ii,j,k,i),j,k,ii
594 write (88,'(f20.10,5h e0(,i1,1h))')
600 write (88,'(a)') restyp(iloctyp(i))
602 write (88,'(f20.15)') b(j,i)
610 & 'Electrostatic interaction constants (asterisk: optimized):'
611 write (iout,'(1x,a,1x,a,7x,a,15x,a,15x,a,13x,a)')
612 & 'IT','JT','EPP','RPP','ELPP6','ELPP3'
615 write(iout,'(2i3,4(1pe15.4,1x,a1,1x))')i,j,
616 & epp(i,j),aster(mask_elec(i,j,1)),
617 & rpp(i,j),aster(mask_elec(i,j,2)),
618 & elpp6(i,j),aster(mask_elec(i,j,3)),
619 & elpp3(i,j),aster(mask_elec(i,j,4))
623 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
624 & 'IT','JT','APP','BPP','AEL6','AEL3'
627 write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
628 & ael6(i,j),ael3(i,j)
631 C Write electrostatic interaction parameters to a file
632 if (nparmset.eq.1) then
633 open(88,file="electr_opt.parm."//prefix(:ilen(prefix)),
636 open(88,file="electr_opt.parm."//prefix(:ilen(prefix))//
640 write(88,'(4f15.10,5x,a)')((epp(i,j),j=1,2),i=1,2),"! EPP"
641 write(88,'(4f15.10,5x,a)')((rpp(i,j),j=1,2),i=1,2),"! RPP"
642 write(88,'(4f15.10,5x,a)')((elpp6(i,j),j=1,2),i=1,2),"! ELPP6"
643 write(88,'(4f15.10,5x,a)')((elpp3(i,j),j=1,2),i=1,2),"! ELPP3"
649 & "Parameters of SC-p interactions (star: optimized):"
650 write (iout,'(t14,a,t34,a)') "TYPE 1","TYPE 2"
651 write (iout,'(8a)') "SC TYPE"," EPS_SCP "," R_SCP ",
652 & " EPS_SCP "," R_SCP "
654 ast11 = aster(mask_scp(i,1,1))
655 ast12 = aster(mask_scp(i,1,2))
656 ast21 = aster(mask_scp(i,2,1))
657 ast22 = aster(mask_scp(i,2,2))
658 if (mask_scp(i,1,1).eq.0) ast11 = aster(mask_scp(0,1,1))
659 if (mask_scp(i,1,2).eq.0) ast12 = aster(mask_scp(0,1,2))
660 if (mask_scp(i,2,1).eq.0) ast21 = aster(mask_scp(0,2,1))
661 if (mask_scp(i,2,2).eq.0) ast22 = aster(mask_scp(0,2,2))
662 write (iout,'(2x,A3,2x,4(f8.3,1x,a1))') restyp(i),
663 & eps_scp(i,1),ast11,rscp(i,1),ast12,eps_scp(i,2),ast21,
666 C Write SCp parameters to a file
667 if (nparmset.eq.1) then
668 open(88,file="scp_opt.parm."//prefix(:ilen(prefix)),
671 open(88,file="scp_opt.parm."//prefix(:ilen(prefix))//
676 write(88,'(4f15.10,5x,a)') eps_scp(i,1),rscp(i,1),
677 & eps_scp(i,2),rscp(i,2),restyp(i)
681 50 format(bz,15(f8.5,1x,a1))
682 60 format(bz,100f10.5)
684 c------------------------------------------------------------------------------
685 subroutine write_zscore
688 include "DIMENSIONS.ZSCOPT"
689 include "COMMON.ENERGIES"
690 include "COMMON.CLASSES"
691 include "COMMON.PROTNAME"
692 include "COMMON.VMCPAR"
693 include "COMMON.IOUNITS"
694 include "COMMON.WEIGHTS"
695 include "COMMON.WEIGHTDER"
696 include "COMMON.COMPAR"
697 include "COMMON.OPTIM"
698 include "COMMON.THERMAL"
699 integer i,j,k,iprot,ib
700 character*6 namnat(5) /"angrms","Q","rmsd","rgy","sign"/
706 write (iout,'(a,2x,a)') "Protein",
707 & protname(iprot)(:ilen(protname(iprot)))
710 write (iout,'(a)')"Maximum likelihood."
712 do ib=1,nbeta(1,iprot)
713 write (iout,'(f8.1,f10.5)')
714 & 1.0d0/(Rgas*betaT(ib,1,iprot)),sumlik(ib,iprot)
719 write (iout,'(a)')"Specific heat."
721 write (iout,'(a,f10.5)') "Baseline",heatbase(iprot)
722 do ib=1,nbeta(2,iprot)
723 write (iout,'(f8.1,3f10.5)')
724 & 1.0d0/(Rgas*betaT(ib,2,iprot)),
725 & zvtot(ib,iprot),target_cv(ib,iprot),weiCv(ib,iprot)
728 if (natlike(iprot).gt.0) write (iout,'(a)')"Nativelikeness."
729 do i=1,natlike(iprot)
730 write (iout,*) "Property",i,numnat(i,iprot)
731 if (natdim(i,iprot).eq.1) then
732 do ib=1,nbeta(i+2,iprot)
733 write (iout,'(f7.2,2f10.5)')
734 & 1.0d0/(Rgas*betaT(ib,i+2,iprot)),nuave(1,ib,i,iprot),
735 & nuexp(1,ib,i,iprot),weinat(1,ib,i,iprot)
738 do ib=1,nbeta(i+2,iprot)
739 write (iout,'("Temperature",f7.2," NATDIM",i5)')
740 & 1.0d0/(Rgas*betaT(ib,i+2,iprot)),natdim(i,iprot)
741 do k=1,natdim(i,iprot)
742 write (iout,'(3f10.5)') nuave(k,ib,i,iprot),
743 & nuexp(k,ib,i,iprot),weinat(k,ib,i,iprot)
753 c------------------------------------------------------------------------------
754 subroutine write_prot(ii,przed)
756 include "DIMENSIONS.ZSCOPT"
757 include "COMMON.ENERGIES"
758 include "COMMON.CLASSES"
759 include "COMMON.PROTNAME"
760 include "COMMON.VMCPAR"
761 include "COMMON.IOUNITS"
762 include "COMMON.COMPAR"
768 nazwa=przed(:ilen(przed))//'.'//prefix(:ilen(prefix))
769 & //'.'//protname(ii)
770 open(unit=istat,file=nazwa)
772 write(istat,'(i10,3f15.3,20e15.5)')
773 & i,e_total(i,ii),eini(i,ii),entfac(i,ii),
774 & (Ptab(i,j,ii),j=1,nbeta(1,ii))