subroutine maxlikopt(nvarr,x) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" #endif include "COMMON.NAMES" include "COMMON.WEIGHTS" include "COMMON.VMCPAR" include "COMMON.OPTIM" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" integer i,j,k,iprot,nf,ii,inn,n,nn,indn(max_ene) logical solved,all_satisfied,not_done double precision ran_number,f,f0,x(max_paropt),x0(max_paropt), & viol integer iran_num,nvarr double precision tcpu,t0,t1,t0w,t1w character*4 liczba integer ilen,lenpre external ilen lenpre=ilen(prefix) #ifdef MPI if (me.eq.master) then if (nparmset.gt.0) then write (liczba,'(bz,i4.4)') myparm restartname=prefix(:lenpre)//"_par"//liczba//'.restart' else restartname=prefix(:lenpre)//'.restart' endif endif #else restartname=prefix(:lenpre)//'.restart' #endif print *,"Start solving inequalities" print *,"Mode",mode if (mode.eq.1) then c Only evaluate initial energies t0=tcpu() #ifdef MPI t0w=MPI_Wtime() #endif call minimize_setup(nvarr,x) t1=tcpu() #ifdef MPI write (iout,*) "CPU time for setup:",t1-t0," sec." t1w=MPI_Wtime() write (iout,*) "Wall clock time for setup:",t1w-t0w," sec." t0w=t1w #endif write (iout,*) do i=1,nvarr x0(i)=x(i) enddo call print_minimize_results(nvarr) return endif if (mode.ge.3) then t0=tcpu() #ifdef MPI t0w=MPI_Wtime() #endif call minimize_setup(nvarr,x) t1=tcpu() #ifdef MPI write (iout,*) "CPU time for setup:",t1-t0," sec." t1w=MPI_Wtime() write (iout,*) "Wall clock time for setup:",t1w-t0w," sec." t0w=t1w #endif t0=t1 print *,"Calling minimize" call minimize(nvarr,x,f0) t1=tcpu() write (iout,*) "CPU time for minimization:",t1-t0," sec." #ifdef MPI t1w=MPI_Wtime() write (iout,*) "Wall clock time for minimization:",t1w-t0w #endif write (iout,*) do i=1,nvarr x0(i)=x(i) enddo c if (nmcm.gt.0) then c write (iout,*) "Start MCM procedure" c call mcm c endif call print_minimize_results(nvarr) return endif c c Scan weight space c if (mode.eq.2) then print *,"Calling scan" call minimize_setup(nvarr,x) call scan(nvarr,x,f0,.true.) do i=1,nvarr xm(i)=x(i) enddo call print_minimize_results(nvarr) endif 50 format(bz,15f8.5) 60 format(bz,100f10.5) return end c----------------------------------------------------------------------------- subroutine minimize_setup(nvar,x) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,TAG,STATUS(MPI_STATUS_SIZE) include "COMMON.MPI" #endif include "COMMON.IOUNITS" include "COMMON.CLASSES" include "COMMON.NAMES" include "COMMON.OPTIM" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.VMCPAR" include "COMMON.ENERGIES" include "COMMON.TIME1" double precision viol integer i,j,k,ii,ib,ibatch,iprot,nvar,nf double precision x(max_paropt),g(max_paropt),rdum,chisquare_force external chisquare_force logical lprn lprn=.true. #ifdef DEBUG write (iout,*) "Entered MINIMIZE_SETUP" call flush(iout) #endif call func1(nvar,3,x_orig) rdum = chisquare_force(nvar,x,g,1) do i=1,n_ene ww_orig(i)=ww(i) enddo call x2w(nvar,x) #ifdef DEBUG write (iout,*) "Initial X",(x(i),i=1,nvar) #endif c write (iout,*) "X_orig",(x_orig(i),i=1,nvar) c write (iout,*) "ww_oorig",(ww_oorig(i),i=1,n_ene) c call flush(iout) call func1(nvar,3,x) rdum = chisquare_force(nvar,x,g,1) #ifdef DEBUG write (iout,*) "x",(x(i),i=1,nvar) #endif write (iout,*) write(iout,'(10(1x,a6,1x))')(wname(print_order(i))(:6),i=1,n_ene) write(iout,40)(ww(print_order(i)),i=1,n_ene) write(iout,*)'-----------------------------------' call write_zscore call flush(iout) return 40 format(10(f7.4,1x)) end c------------------------------------------------------------------------------ subroutine print_minimize_results(nvarr) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer ierror,status(MPI_STATUS_SIZE),tag include "COMMON.MPI" #endif include "COMMON.CONTROL" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.PROTNAME" include "COMMON.ENERGIES" include "COMMON.CLASSES" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.OPTIM" include "COMMON.TIME1" integer nf,nvarr,i,ib,ic,ii,num_nat_tot,iprot,igather double precision x(max_ene),g(max_ene) double precision viol,fdum,rdum,chisquare_force external chisquare_force character*4 liczba logical lprn lprn=.true. c write (iout,*) "print_minimize_results" c call flush(iout) write(iout,*) '----------------------' cutoffeval=.false. call x2w(nvarr,xm(1)) c write (iout,*) "Calling func1" c call flush(iout) call func1(nvarr,2,xm(1)) rdum = chisquare_force(nvarr,x,g,3) c write (iout,*) "Calling write_params" c call flush(iout) call write_params(nvarr,nf,xm(1)) c write(iout,*) "After write_params" c call flush(iout) if (out_newe) then do ii=1,nprot if (maxlik(ii)) then if (nparmset.eq.1) then call write_prot(ii,'NewE_all') else write (liczba,'(bz,i4.4)') myparm call write_prot(ii,'NewE_all'//liczba) endif c do ic=1,nclass(ii) c call make_distrib(iclass(ic,ii),ii) c enddo endif enddo endif if (out_forces) then do ii=1,nprot if (fmatch(ii)) call write_forces(ii) enddo endif return end c------------------------------------------------------------------------------ subroutine write_params(nvar,nf,x) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.CONTROL" include "COMMON.WEIGHTS" include "COMMON.NAMES" include "COMMON.INTERACT" include "COMMON.TORSION" include "COMMON.LOCAL" include "COMMON.SCCOR" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.FFIELD" include "COMMON.OPTIM" integer nf,nvar double precision x(max_paropt) integer ilen external ilen integer i,ii,j,l,ll,jj,pj,jstart,jend,it1,it2,k,iblock integer itor2typ(-2:2) /-20,9,10,9,20/ character*1 aster(0:1) /" ","*"/,ast11,ast12,ast21,ast22 character*4 liczba character*1 toronelet(-2:2) /"p","a","G","A","P"/ logical lprint call x2w(nvar,x) call write_zscore write (iout,*) "Weights (asterisk: optimized)" write(iout,'(1x,15(a8,2x))') (wname(print_order(i))(:8),i=1,n_ene) write(iout,50)(ww(print_order(i)), & aster(imask(print_order(i))),i=1,n_ene) write(iout,*) C Write weights in UNRES format for the input file if (nparmset.eq.1) then open(88,file="weights_opt."//prefix(:ilen(prefix)), & status="unknown") else write (liczba,'(bz,i4.4)') myparm open(88,file="weights_opt."//prefix(:ilen(prefix))//"par_" & //liczba, & status="unknown") endif c do i=1,n_ene,5 do i=1,n_ene,4 jstart=i c jend=min0(i+4,n_ene) jend=min0(i+3,n_ene) jj=0 do j=jstart,jend pj=print_order(j) c write(88,'(bz,2a,f7.5," ",$)') c & wname(pj)(:ilen(wname(pj))),'=',ww(pj) c jj=jj+ilen(wname(pj))+9 write(88,'(bz,2a,e11.5," ",$)') & wname(pj)(:ilen(wname(pj))),'=',ww(pj) jj=jj+ilen(wname(pj))+13 enddo write (88,'(a,$)') (" ",j=1,79-jj) write (88,'("&")') enddo write (88,'(bz,2(2a,f7.5," "))') & "CUTOFF","=",7.0d0,"WCORR4","=",0.0d0 close(88) if (mod_side .or. mod_side_other) then if (nparmset.eq.1) then open(88, & file="sc_"//POT(:ilen(POT))//"_opt."//prefix(:ilen(prefix)), & status="unknown") else open(88, & file="sc_"//POT(:ilen(POT))//"_opt."//prefix(:ilen(prefix))// & "_par"//liczba, & status="unknown") endif lprint=.true. write(88,'(2i5)') ipot,expon if (ipot.gt.2) then do i=1,ntyp chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0) enddo endif write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot), & ', exponents are ',expon,2*expon goto (10,20,30,30,40) ipot C----------------------- LJ potential --------------------------------- 10 do i=1,ntyp write (88,'(4f20.10)')(eps(i,j),j=i,ntyp) write (iout,*) enddo write (88,'(4f20.10)') (sigma0(i),i=1,ntyp) write (88,*) if (lprint) then write (iout,'(/a/)') 'Parameters of the LJ potential:' write (iout,'(a/)') 'The epsilon array:' call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) write (iout,'(/a)') 'One-body parameters:' write (iout,'(a,4x,a)') 'residue','sigma' write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp) endif goto 55 C----------------------- LJK potential -------------------------------- 20 do i=1,ntyp write (88,'(4f20.10)')(eps(i,j),j=i,ntyp) write (iout,*) enddo write (88,'(4f20.10)') (sigma0(i),i=1,ntyp) write (88,*) write (88,'(4f20.10)') (rr0(i),i=1,ntyp) write (88,*) if (lprint) then write (iout,'(/a/)') 'Parameters of the LJK potential:' write (iout,'(a/)') 'The epsilon array:' call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) write (iout,'(/a)') 'One-body parameters:' write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 ' write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i), & i=1,ntyp) endif goto 55 C---------------------- GB or BP potential ----------------------------- 30 do i=1,ntyp write (88,'(4f20.15)') (eps(i,j),j=i,ntyp) write (88,*) enddo write (88,'(4f20.15)')(sigma0(i),i=1,ntyp) write (88,*) write (88,'(4f20.15)')(sigii(i),i=1,ntyp) write (88,*) write (88,'(4f20.15)')(chip0(i),i=1,ntyp) write (88,*) write (88,'(4f20.15)')(alp(i),i=1,ntyp) write (88,*) C For the GB potential convert sigma'**2 into chi' if (lprint) then if (ipot.eq.3) then write (iout,'(/a/)') 'Parameters of the BP potential:' else write (iout,'(/a/)') 'Parameters of the GB potential:' endif write (iout,'(a/)') 'The epsilon array:' call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) write (iout,'(/a)') 'One-body parameters (asterisk: optimized):' write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2', & ' chip ',' alph ' do i=1,ntyp write (iout,'(a3,6x,f10.5,1x,a1,3(f8.5,1x,a1))') restyp(i), & sigma0(i),aster(mask_sigma(1,i)),sigii(i),aster(mask_sigma(2,i)), & chip(i),aster(mask_sigma(3,i)),alp(i),aster(mask_sigma(4,i)) enddo endif goto 55 C--------------------- GBV potential ----------------------------------- 40 do i=1,ntyp write (88,'(4f20.15)') (eps(i,j),j=i,ntyp) write (88,*) enddo write (88,'(4f20.15)')(sigma0(i),i=1,ntyp) write (88,*) write (88,'(4f20.15)')(rr0(i),i=1,ntyp) write (88,*) write (88,'(4f20.15)')(sigii(i),i=1,ntyp) write (88,*) write (88,'(4f20.15)')(chip0(i),i=1,ntyp) write (88,*) write (88,'(4f20.15)')(alp(i),i=1,ntyp) write (88,*) if (lprint) then write (iout,'(/a/)') 'Parameters of the GBV potential:' write (iout,'(a/)') 'The epsilon array:' call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) write (iout,'(/a)') 'One-body parameters:' write (iout,'(a,4x,6a)') 'residue',' sigma ',' r0 ', & 's||/s_|_^2',' chip0 ',' chip ',' alph ' do i=1,ntyp write (iout,'(a3,6x,f10.5,1x,a1,5(f8.5,1x,a1))') restyp(i), & sigma0(i),aster(mask_sigma(1,i)),rr0(i),aster(mask_sigma(5,i)), & sigii(i),aster(mask_sigma(2,i)),chip0(i),aster(mask_sigma(3,i)), & chip(i),aster(mask_sigma(3,i)),alp(i),aster(mask_sigma(4,i)) enddo endif 55 continue C----------------------------------------------------------------------- close(88) endif #ifdef NEWCORR if (mod_tor .or. tor_mode.eq.2) then #else if (mod_tor) then #endif write (iout,'(a)') "Optimized weights of the torsional terms." do iblock=1,2 do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 if (mask_tor(0,j,i,iblock).gt.0) then write (iout,'(i4,2a4,f10.5)') 0,restyp(iloctyp(i)), & restyp(iloctyp(j)),weitor(0,j,i,iblock) endif enddo enddo enddo if (nparmset.eq.1) then open(88,file="tor_opt.parm."//prefix(:ilen(prefix)), & status="unknown") write (88,'(i1,7h *** ,2a,7h *** )') ntortyp, & "Optimized torsional parameters ",prefix(:ilen(prefix)) write (88,'(40i2)') (itortyp(i),i=1,ntyp) if (tor_mode.eq.0) then write (iout,'(/a/)') 'Torsional constants:' do iblock=1,2 do i=0,ntortyp-1 do j=0,ntortyp-1 write (iout,*) 'ityp',i,' jtyp',j write (iout,*) 'Fourier constants' do k=1,nterm(i,j,iblock) write (iout,'(2(1pe15.5))') v1(k,i,j,iblock), & v2(k,i,j,iblock) enddo write (iout,*) 'Lorenz constants' do k=1,nlor(i,j,iblock) write (iout,'(3(1pe15.5))') & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) enddo enddo enddo enddo do iblock=1,2 do i=0,ntortyp-1 do j=-ntortyp+1,ntortyp-1 write (88,'(2i2,16h ************ ,a1,1h-,a1)') & nterm(i,j,iblock),nlor(i,j,iblock),onelet(iloctyp(i)), & onelet(iloctyp(j)) do k=1,nterm(i,j,iblock) write (88,'(i6,13x,2(1pe18.5))') k, & v1(k,i,j,iblock)*weitor(0,i,j,iblock),v2(k,i,j, & iblock)*weitor(0,i,j,iblock) enddo do k=1,nlor(i,j,iblock) write (88,'(i6,13x,3(1pe18.5))') & k,vlor1(k,i,j)*weitor(k,i,j,iblock), & vlor2(k,i,j)*weitor(k,i,j,iblock), & vlor3(k,i,j)*weitor(k,i,j,iblock) enddo enddo enddo enddo else c Print valence-torsional parameters write (iout,'(a)') & "Parameters of the valence-torsional potentials" do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 write (iout,'(3a)')"Type ",onelet(iloctyp(i)),onelet(iloctyp(j)) write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc" do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) write (iout,'(3i5,2f15.4)') & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i) enddo enddo enddo enddo enddo do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 write (88,'(16h ************ ,a1,1h-,a1)') & onelet(iloctyp(i)),onelet(iloctyp(j)) write (88,'(2i5)') nterm_kcc(j,i),nterm_kcc_Tb(j,i) do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) write (88,'(3i5,2(1pe15.5))') k,l-1,ll-1, & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i) enddo enddo enddo enddo enddo close(88) endif endif endif if (mod_sccor) then write (iout,'(a)') "Optimized weights of the sccor terms." do i=1,nsccortyp do j=1,nsccortyp do l=1,3 if (mask_tor(l,j,i,iblock).gt.0) then write (iout,'(i4,2a4,f10.5)') l,restyp(isccortyp(j)), & restyp(isccortyp(i)),weitor(l,j,i,1) endif enddo enddo enddo if (nparmset.eq.1) then open(88,file="sccor_opt.parm."//prefix(:ilen(prefix)), & status="unknown") write (88,'(i2,7h *** ,2a,7h *** )') nsccortyp, & "Optimized local correlation parameters ",prefix(:ilen(prefix)) write (88,'(20i4)') (isccortyp(i),i=1,ntyp) do l=1,3 do i=1,nsccortyp do j=1,nsccortyp write (88,'(2i4)') nterm_sccor(i,j),nlor_sccor(i,j) write (88,'(12(1h*),1x,a1,1h-,a1)') onelet(i),onelet(j) do k=1,nterm_sccor(i,j) write (88,'(i6,13x,2(1pe18.5))') k, & v1sccor(k,l,i,j)*weitor(l,i,j,1), & v2sccor(k,l,i,j)*weitor(l,i,j,1) enddo enddo enddo enddo close(88) endif endif C----------------------------------------------------------------------- c 7/8/17 AL: Optimization of the bending parameters if (mod_ang) then write(iout,*) & "Optimized angle potential coefficients (same for D-types)" do i=0,nthetyp-1 if (mask_ang(i).eq.0) cycle write (iout,*) "Type: ",onelet(iloctyp(i)) do j=1,nbend_kcc_TB(i) write (iout,'(i5,f15.5)') j,v1bend_chyb(j,i) enddo enddo open(88,file="theta_opt.parm."//prefix(:ilen(prefix)), & status="unknown") write (88,'(i2)') nthetyp do i=-nthetyp+1,nthetyp-1 write (88,'(i2,1x,12(1h*),1x,a1)') & nbend_kcc_TB(i),onelet(iloctyp(i)) do j=0,nbend_kcc_TB(i) write (88,'(i5,f20.10)') j,v1bend_chyb(j,i) enddo enddo close(88) endif C----------------------------------------------------------------------- if (mod_fourier(nloctyp).gt.0) then #ifdef NEWCORR write (iout,'(2a)') & "Fourier coefficients of new (angle-dependent) cumulants", & " (asterisk: optimized)" do i=0,nloctyp-1 write (iout,*) "Type: ",onelet(iloctyp(i)) write (iout,*) "Coefficients of the expansion of B1" do j=1,2 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of B2" do j=1,2 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of C" write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3) write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of D" write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3) write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of E" write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3) do j=1,2 do k=1,2 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2) enddo enddo enddo IF (SPLIT_FOURIERTOR) THEN write (iout,'(2a)') & "Fourier coefficients of new (angle-dependent) cumulants", & " for torsional potentials (asterisk: optimized)" do i=0,nloctyp-1 write (iout,*) "Type: ",onelet(iloctyp(i)) write (iout,*) "Coefficients of the expansion of B1" do j=1,2 write (iout,'(3hB1(,i1,1h),3f10.5)') & j,(bnew1tor(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of B2" do j=1,2 write (iout,'(3hB2(,i1,1h),3f10.5)') & j,(bnew2tor(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of C" write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3) write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of D" write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3) write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of E" write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3) do j=1,2 do k=1,2 write (iout,'(1hE,2i1,2f10.5)') & j,k,(eenewtor(l,j,k,i),l=1,2) enddo enddo enddo ENDIF #else write (iout,'(2a)') & "Fourier coefficients of old angle-independent cumulants", & " (asterisk: optimized)" do i=0,nloctyp-1 write (iout,*) 'Type ',restyp(iloctyp(i)) write (iout,'(a,i2,a,f10.5,1x,a)') & ('b(',j,')=',b(j,i),aster(mask_fourier(j,i)),j=1,13) enddo #endif c Write a file with parameters for UNRES if (nparmset.eq.1) then open(88,file="fourier_opt.parm."//prefix(:ilen(prefix)), & status="unknown") else open(88,file="fourier_opt.parm."//prefix(:ilen(prefix))// & "_par"//liczba, & status="unknown") endif #ifdef NEWCORR IF (SPLIT_FOURIERTOR) THEN write (88,'(i5,5x,a)') -nloctyp, & "# Number of local interaction types" ELSE write (88,'(i5,5x,a)') nloctyp, & "# Number of local interaction types" ENDIF write (88,'(40i2)') (itype2loc(i),i=1,ntyp) write (88,'(20i4)') (iloctyp(i),i=0,nloctyp) do i=0,nloctyp-1 write (88,'(a)') restyp(iloctyp(i)) do ii=1,3 do j=1,2 write (88,'(f20.10,4h b1,i1,1h(,i1,1h))') & bnew1(ii,j,i),j,ii enddo enddo do ii=1,3 do j=1,2 write (88,'(f20.10,4h b2,i1,1h(,i1,1h))') & bnew2(ii,j,i),j,ii enddo enddo do ii=1,3 do j=1,2 write (88,'(f20.10,4h c1,i1,1h(,i1,1h))') & 2*ccnew(ii,j,i),j,ii enddo enddo do ii=1,3 do j=1,2 write (88,'(f20.10,4h d1,i1,1h(,i1,1h))') & 2*ddnew(ii,j,i),j,ii enddo enddo do ii=1,2 do j=1,2 do k=1,2 write (88,'(f20.10,3h e,2i1,1h(,i1,1h))') & eenew(ii,j,k,i),j,k,ii enddo enddo enddo do ii=1,3 write (88,'(f20.10,5h e0(,i1,1h))') & e0new(ii,i),ii enddo enddo IF (SPLIT_FOURIERTOR) THEN do i=0,nloctyp-1 write (88,'(a)') restyp(iloctyp(i)) do ii=1,3 do j=1,2 write (88,'(f20.10,4h b1,i1,1h(,i1,1h))') & bnew1tor(ii,j,i),j,ii enddo enddo do ii=1,3 do j=1,2 write (88,'(f20.10,4h b2,i1,1h(,i1,1h))') & bnew2tor(ii,j,i),j,ii enddo enddo do ii=1,3 do j=1,2 write (88,'(f20.10,4h c1,i1,1h(,i1,1h))') & 2*ccnewtor(ii,j,i),j,ii enddo enddo do ii=1,3 do j=1,2 write (88,'(f20.10,4h d1,i1,1h(,i1,1h))') & 2*ddnewtor(ii,j,i),j,ii enddo enddo do ii=1,2 do j=1,2 do k=1,2 write (88,'(f20.10,3h e,2i1,1h(,i1,1h))') & eenewtor(ii,j,k,i),j,k,ii enddo enddo enddo do ii=1,3 write (88,'(f20.10,5h e0(,i1,1h))') & e0newtor(ii,i),ii enddo enddo ENDIF #else write (88,'(i5,5x,a)') nloctyp, & "# Number of local interaction types" do i=1,nloctyp write (88,'(a)') restyp(iloctyp(i)) do j=1,13 write (88,'(f20.15)') b(j,i) enddo enddo #endif close(88) endif if (mod_elec) then write (iout,'(/a)') & 'Electrostatic interaction constants (asterisk: optimized):' write (iout,'(1x,a,1x,a,7x,a,15x,a,15x,a,13x,a)') & 'IT','JT','EPP','RPP','ELPP6','ELPP3' do i=1,2 do j=1,2 write(iout,'(2i3,4(1pe15.4,1x,a1,1x))')i,j, & epp(i,j),aster(mask_elec(i,j,1)), & rpp(i,j),aster(mask_elec(i,j,2)), & elpp6(i,j),aster(mask_elec(i,j,3)), & elpp3(i,j),aster(mask_elec(i,j,4)) enddo enddo write (iout,*) write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') & 'IT','JT','APP','BPP','AEL6','AEL3' do i=1,2 do j=1,2 write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j), & ael6(i,j),ael3(i,j) enddo enddo C Write electrostatic interaction parameters to a file if (nparmset.eq.1) then open(88,file="electr_opt.parm."//prefix(:ilen(prefix)), & status="unknown") else open(88,file="electr_opt.parm."//prefix(:ilen(prefix))// & "_par"//liczba, & status="unknown") endif write(88,'(4f15.10,5x,a)')((epp(i,j),j=1,2),i=1,2),"! EPP" write(88,'(4f15.10,5x,a)')((rpp(i,j),j=1,2),i=1,2),"! RPP" write(88,'(4f15.10,5x,a)')((elpp6(i,j),j=1,2),i=1,2),"! ELPP6" write(88,'(4f15.10,5x,a)')((elpp3(i,j),j=1,2),i=1,2),"! ELPP3" close(88) endif if (mod_scp) then write (iout,*) write (iout,*) & "Parameters of SC-p interactions (star: optimized):" write (iout,'(t14,a,t34,a)') "TYPE 1","TYPE 2" write (iout,'(8a)') "SC TYPE"," EPS_SCP "," R_SCP ", & " EPS_SCP "," R_SCP " do i=1,20 ast11 = aster(mask_scp(i,1,1)) ast12 = aster(mask_scp(i,1,2)) ast21 = aster(mask_scp(i,2,1)) ast22 = aster(mask_scp(i,2,2)) if (mask_scp(i,1,1).eq.0) ast11 = aster(mask_scp(0,1,1)) if (mask_scp(i,1,2).eq.0) ast12 = aster(mask_scp(0,1,2)) if (mask_scp(i,2,1).eq.0) ast21 = aster(mask_scp(0,2,1)) if (mask_scp(i,2,2).eq.0) ast22 = aster(mask_scp(0,2,2)) write (iout,'(2x,A3,2x,4(f8.3,1x,a1))') restyp(i), & eps_scp(i,1),ast11,rscp(i,1),ast12,eps_scp(i,2),ast21, & rscp(i,2),ast22 enddo C Write SCp parameters to a file if (nparmset.eq.1) then open(88,file="scp_opt.parm."//prefix(:ilen(prefix)), & status="unknown") else open(88,file="scp_opt.parm."//prefix(:ilen(prefix))// & "_par"//liczba, & status="unknown") endif do i=1,ntyp write(88,'(4f15.10,5x,a)') eps_scp(i,1),rscp(i,1), & eps_scp(i,2),rscp(i,2),restyp(i) enddo close(88) endif 50 format(bz,15(f8.5,1x,a1)) 60 format(bz,100f10.5) end c------------------------------------------------------------------------------ subroutine write_zscore implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.ENERGIES" include "COMMON.CLASSES" include "COMMON.PROTNAME" include "COMMON.VMCPAR" include "COMMON.IOUNITS" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.COMPAR" include "COMMON.OPTIM" include "COMMON.THERMAL" include "COMMON.FMATCH" integer i,j,k,iprot,ib character*6 namnat(5) /"angrms","Q","rmsd","rgy","sign"/ integer ilen external ilen do iprot=1,nprot write (iout,*) write (iout,'(a,2x,a)') "Protein", & protname(iprot)(:ilen(protname(iprot))) if (maxlik(iprot)) then write (iout,*) write (iout,'(a)')"Maximum likelihood." write (iout,*) do ib=1,nbeta(1,iprot) write (iout,'(f8.1,f10.5)') & 1.0d0/(Rgas*betaT(ib,1,iprot)),sumlik(ib,iprot) enddo write (iout,*) endif if (nbeta(2,iprot).gt.0) then write (iout,*) write (iout,'(a)')"Specific heat." write (iout,*) write (iout,'(a,f10.5)') "Baseline",heatbase(iprot) do ib=1,nbeta(2,iprot) write (iout,'(f8.1,3f10.5)') & 1.0d0/(Rgas*betaT(ib,2,iprot)), & zvtot(ib,iprot),target_cv(ib,iprot),weiCv(ib,iprot) enddo endif if (natlike(iprot).gt.0) write (iout,'(a)')"Nativelikeness." do i=1,natlike(iprot) write (iout,*) "Property",i,numnat(i,iprot) if (natdim(i,iprot).eq.1) then do ib=1,nbeta(i+2,iprot) write (iout,'(f7.2,2f10.5)') & 1.0d0/(Rgas*betaT(ib,i+2,iprot)),nuave(1,ib,i,iprot), & nuexp(1,ib,i,iprot),weinat(1,ib,i,iprot) enddo else do ib=1,nbeta(i+2,iprot) write (iout,'("Temperature",f7.2," NATDIM",i5)') & 1.0d0/(Rgas*betaT(ib,i+2,iprot)),natdim(i,iprot) do k=1,natdim(i,iprot) write (iout,'(3f10.5)') nuave(k,ib,i,iprot), & nuexp(k,ib,i,iprot),weinat(k,ib,i,iprot) enddo enddo endif enddo if (fmatch(iprot)) then write (iout,'(a32,f10.5,a)') & "Force mean deviation:",dsqrt(chisquare(iprot))," kcal/(mol*A)" write (iout,'(a32,f10.5,a)') & "MD force deviation:",dsqrt(varforce(iprot))," kcal/(mol*A)" write (iout,'(a32,f10.5,a)') & "Correlation coeff.:", & rcorr(iprot) endif enddo write (iout,*) return end c------------------------------------------------------------------------------ subroutine write_forces(iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR include "COMMON.MPI" #endif include "COMMON.ENERGIES" include "COMMON.CLASSES" include "COMMON.PROTNAME" include "COMMON.VMCPAR" include "COMMON.IOUNITS" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.COMPAR" include "COMMON.OPTIM" include "COMMON.THERMAL" include "COMMON.FMATCH" include "COMMON.NAMES" include "COMMON.INTERACT" include "COMMON.CHAIN" character*128 fnam integer i,j,k,jj,iprot integer ilen external ilen character*4 liczba if (nparmset.gt.1) then write (liczba,'(bz,i4.4)') myparm fnam = prefix(:ilen(prefix))// & protname(iprot)(:ilen(protname(iprot)))//liczba//"_forces.plt" else fnam = prefix(:ilen(prefix))// & protname(iprot)(:ilen(protname(iprot)))//"_forces.plt" endif open (istat,file=fnam,status="unknown") write (iout,*) "write_forces iprot ntot",iprot,ntot_work_MD(iprot) do i=1,ntot_work_MD(iprot) do j=nnt,nct write(istat,'(1h"a1,i2.2,1h",i3,i6,3f10.4,5x,3f10.4,5x,3f10.4)') & onelet(itype(j)),j-1,j,i,(CGforce_MD(k,j,i,iprot),k=1,3), & (force(k,j,i,iprot),k=1,3),(CGforce_MD_ave(k,j,i,iprot),k=1,3) c write(iout,'(1h"a1,i2.2,1h",i3,i6,3f10.4,5x,3f10.4)') c & onelet(itype(j)),j-1,j,i,(CGforce_MD(k,j,i,iprot),k=1,3), c & (force(k,j,i,iprot),k=1,3) enddo jj=nres do j=nnt,nct if (itype(j).ne.10. and. itype(j).ne.ntyp1) then jj=jj+1 write(istat,'(1h"a1,i2.2,1h",i3,i6,3f10.4,5x,3f10.4,5x,3f10.4)') & onelet(itype(j)),j-1,jj,i,(CGforce_MD(k,jj,i,iprot),k=1,3), & (force(k,jj,i,iprot),k=1,3),(CGforce_MD_ave(k,jj,i,iprot),k=1,3) c write(iout,'(1h"a1,i2.2,1h",i3,i6,3f10.4,5x,3f10.4)') c & onelet(itype(j)),j-1,jj,i,(CGforce_MD(k,jj,i,iprot), c & k=1,3),(force(k,jj,i,iprot),k=1,3) endif enddo enddo close(istat) return end c------------------------------------------------------------------------------ subroutine write_prot(ii,przed) implicit none include "DIMENSIONS.ZSCOPT" include "COMMON.ENERGIES" include "COMMON.CLASSES" include "COMMON.PROTNAME" include "COMMON.VMCPAR" include "COMMON.IOUNITS" include "COMMON.COMPAR" integer i,ii,j integer ilen external ilen character*(*) przed character*64 nazwa nazwa=przed(:ilen(przed))//'.'//prefix(:ilen(prefix)) & //'.'//protname(ii) open(unit=istat,file=nazwa) do i=1,ntot_work(ii) write(istat,'(i10,3f15.3,20e15.5)') & i,e_total(i,ii),eini(i,ii),entfac(i,ii), & (Ptab(i,j,ii),j=1,nbeta(1,ii)) enddo close(istat) return end