source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-PDB/maxlikopt.F

   1       subroutine maxlikopt(nvarr,x)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5 #ifdef MPI
   6       include "mpif.h"
   7       include "COMMON.MPI"
   8 #endif
   9       include "COMMON.NAMES"
  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),x0(max_paropt),
  19      &  viol
  20       integer iran_num,nvarr
  21       double precision tcpu,t0,t1,t0w,t1w
  22       character*4 liczba
  23       integer ilen,lenpre
  24       external ilen
  25
  26       lenpre=ilen(prefix)
  27 #ifdef MPI
  28       if (me.eq.master) then
  29         if (nparmset.gt.0) then
  30           write (liczba,'(bz,i4.4)') myparm
  31           restartname=prefix(:lenpre)//"_par"//liczba//'.restart'
  32         else
  33           restartname=prefix(:lenpre)//'.restart'
  34         endif
  35       endif
  36 #else
  37       restartname=prefix(:lenpre)//'.restart'
  38 #endif
  39       print *,"Start solving inequalities"
  40       print *,"Mode",mode
  41       if (mode.eq.1) then
  42 c Only evaluate initial energies
  43         t0=tcpu()
  44 #ifdef MPI
  45         t0w=MPI_Wtime()
  46 #endif
  47         call minimize_setup(nvarr,x)
  48         t1=tcpu()
  49 #ifdef MPI
  50         write (iout,*) "CPU time for setup:",t1-t0," sec."
  51         t1w=MPI_Wtime()
  52         write (iout,*) "Wall clock time for setup:",t1w-t0w," sec."
  53         t0w=t1w
  54 #endif
  55         write (iout,*)
  56         do i=1,nvarr
  57           x0(i)=x(i)
  58         enddo
  59         call print_minimize_results(nvarr)
  60         return
  61       endif
  62       if (mode.ge.3) then
  63         t0=tcpu()
  64 #ifdef MPI
  65         t0w=MPI_Wtime()
  66 #endif
  67         call minimize_setup(nvarr,x)
  68         t1=tcpu()
  69 #ifdef MPI
  70         write (iout,*) "CPU time for setup:",t1-t0," sec."
  71         t1w=MPI_Wtime()
  72         write (iout,*) "Wall clock time for setup:",t1w-t0w," sec."
  73         t0w=t1w
  74 #endif
  75         t0=t1
  76         print *,"Calling minimize"
  77         call minimize(nvarr,x,f0)
  78         t1=tcpu()
  79         write (iout,*) "CPU time for minimization:",t1-t0," sec."
  80 #ifdef MPI
  81         t1w=MPI_Wtime()
  82         write (iout,*) "Wall clock time for minimization:",t1w-t0w
  83 #endif
  84         write (iout,*)
  85         do i=1,nvarr
  86           x0(i)=x(i)
  87         enddo
  88 c        if (nmcm.gt.0) then
  89 c          write (iout,*) "Start MCM procedure"
  90 c          call mcm
  91 c        endif
  92         call print_minimize_results(nvarr)
  93         return
  94       endif
  95 c
  96 c Scan weight space
  97 c
  98       if (mode.eq.2) then
  99         print *,"Calling scan"
 100         call minimize_setup(nvarr,x)
 101         call scan(nvarr,x,f0,.true.)
 102         do i=1,nvarr
 103           xm(i)=x(i)
 104         enddo
 105         call print_minimize_results(nvarr)
 106       endif
 107 50    format(bz,15f8.5)
 108 60    format(bz,100f10.5)
 109       return
 110       end
 111 c-----------------------------------------------------------------------------
 112       subroutine minimize_setup(nvar,x)
 113       implicit none
 114       include "DIMENSIONS"
 115       include "DIMENSIONS.ZSCOPT"
 116 #ifdef MPI
 117       include "mpif.h"
 118       integer IERROR,TAG,STATUS(MPI_STATUS_SIZE)
 119       include "COMMON.MPI"
 120 #endif
 121       include "COMMON.IOUNITS"
 122       include "COMMON.CLASSES"
 123       include "COMMON.NAMES"
 124       include "COMMON.OPTIM"
 125       include "COMMON.TORSION"
 126       include "COMMON.WEIGHTS"
 127       include "COMMON.WEIGHTDER"
 128       include "COMMON.VMCPAR"
 129       include "COMMON.ENERGIES"
 130       include "COMMON.TIME1"
 131       double precision g(nvar)
 132       double precision viol
 133       integer i,j,k,ii,ib,ibatch,iprot,nvar,nf
 134       double precision x(max_paropt)
 135       double precision rdif,maxlik_pdb
 136       logical lprn
 137       lprn=.true.
 138 #ifdef DEBUG
 139       write (iout,*) "Entered MINIMIZE_SETUP"
 140       call flush(iout)
 141 #endif
 142       call func1(nvar,3,x_orig)
 143       do i=1,n_ene
 144         ww_orig(i)=ww(i)
 145       enddo
 146       call x2w(nvar,x)
 147 c      write (iout,*) "Initial X",(x(i),i=1,nvar)
 148 c      write (iout,*) "X_orig",(x_orig(i),i=1,nvar)
 149 c      write (iout,*) "ww_oorig",(ww_oorig(i),i=1,n_ene)
 150 c      call flush(iout)
 151       if (torsion_pmf.or.turn_pmf.or.eello_pmf.or.angle_PMF.or.pdbpmf)
 152      &then
 153         call func1(nvar,1,x)
 154         if (torsion_pmf.or.turn_pmf.or.eello_pmf.or.angle_PMF)
 155      &    fpmf = rdif(nvar,x,g,1)
 156         if (pdbpmf) fpdb = maxlik_pdb(nvar,x,g,1)
 157       endif
 158       call func1(nvar,3,x)
 159 #ifdef DEBUG
 160       write (iout,*) "x",(x(i),i=1,nvar)
 161 #endif
 162       write (iout,*)
 163       write(iout,'(10(1x,a6,1x))')(wname(print_order(i))(:6),i=1,n_ene)
 164       write(iout,40)(ww(print_order(i)),i=1,n_ene)
 165       write(iout,*)'-----------------------------------'
 166       call write_zscore
 167       call flush(iout)
 168       return
 169 40    format(10(f7.4,1x))
 170       end
 171 c------------------------------------------------------------------------------
 172       subroutine print_minimize_results(nvarr)
 173       implicit none
 174       include "DIMENSIONS"
 175       include "DIMENSIONS.ZSCOPT"
 176 #ifdef MPI
 177       include "mpif.h"
 178       integer ierror,status(MPI_STATUS_SIZE),tag
 179       include "COMMON.MPI"
 180 #endif
 181       include "COMMON.CONTROL"
 182       include "COMMON.WEIGHTS"
 183       include "COMMON.WEIGHTDER"
 184       include "COMMON.PROTNAME"
 185       include "COMMON.ENERGIES"
 186       include "COMMON.CLASSES"
 187       include "COMMON.IOUNITS"
 188       include "COMMON.VMCPAR"
 189       include "COMMON.OPTIM"
 190       include "COMMON.TIME1"
 191       integer nf,nvarr,i,ib,ic,ii,num_nat_tot,iprot,igather
 192       double precision x(max_ene)
 193       double precision viol
 194       character*4 liczba
 195       logical lprn
 196       lprn=.true.
 197 c      write (iout,*) "print_minimize_results"
 198       call flush(iout)
 199       write(iout,*) '----------------------'
 200       cutoffeval=.false.
 201       call x2w(nvarr,xm(1))
 202 c      write (iout,*) "Calling func1"
 203       call flush(iout)
 204       call func1(nvarr,2,xm(1))
 205 c      write (iout,*) "Calling write_params"
 206       call flush(iout)
 207       call write_params(nvarr,nf,xm(1))
 208 c      write(iout,*) "After write_params"
 209       call flush(iout)
 210       if (out_newe) then
 211       do ii=1,nprot
 212         if (nparmset.eq.1) then
 213         call write_prot(ii,'NewE_all')
 214         else
 215         write (liczba,'(bz,i4.4)') myparm
 216         call write_prot(ii,'NewE_all'//liczba)
 217         endif
 218 c        do ic=1,nclass(ii)
 219 c          call make_distrib(iclass(ic,ii),ii)
 220 c        enddo
 221       enddo
 222       endif
 223
 224       return
 225
 226       end
 227 c------------------------------------------------------------------------------
 228       subroutine write_params(nvar,nf,x)
 229       implicit none
 230       include "DIMENSIONS"
 231       include "DIMENSIONS.ZSCOPT"
 232       include "COMMON.CONTROL"
 233       include "COMMON.WEIGHTS"
 234       include "COMMON.NAMES"
 235       include "COMMON.INTERACT"
 236       include "COMMON.TORSION"
 237       include "COMMON.LOCAL"
 238       include "COMMON.SCCOR"
 239       include "COMMON.CLASSES"
 240       include "COMMON.ENERGIES"
 241       include "COMMON.IOUNITS"
 242       include "COMMON.FFIELD"
 243       include "COMMON.OPTIM"
 244       integer nf,nvar
 245       double precision x(max_paropt)
 246       integer ilen
 247       external ilen
 248       integer i,ii,j,l,ll,jj,pj,jstart,jend,it1,it2,k,iblock
 249       integer itor2typ(-2:2) /-20,9,10,9,20/
 250       character*1 aster(0:1) /" ","*"/,ast11,ast12,ast21,ast22
 251       character*4 liczba
 252       character*1 toronelet(-2:2) /"p","a","G","A","P"/
 253       logical lprint
 254       integer mysign
 255       call x2w(nvar,x)
 256       call write_zscore
 257       write (iout,*) "Weights (asterisk: optimized)"
 258       write(iout,'(1x,15(a8,2x))') (wname(print_order(i))(:8),i=1,n_ene)
 259       write(iout,50)(ww(print_order(i)),
 260      &  aster(mysign(imask(print_order(i)))),i=1,n_ene)
 261       write(iout,*)
 262 C Write weights in UNRES format for the input file
 263       if (nparmset.eq.1) then
 264         open(88,file="weights_opt."//prefix(:ilen(prefix)),
 265      &   status="unknown")
 266       else
 267         write (liczba,'(bz,i4.4)') myparm
 268         open(88,file="weights_opt."//prefix(:ilen(prefix))//"par_"
 269      &   //liczba,
 270      &   status="unknown")
 271       endif
 272 c      do i=1,n_ene,5
 273       do i=1,n_ene,4
 274         jstart=i
 275 c        jend=min0(i+4,n_ene)
 276         jend=min0(i+3,n_ene)
 277         jj=0
 278         do j=jstart,jend
 279           pj=print_order(j)
 280 c          write(88,'(bz,2a,f7.5," ",$)')
 281 c     &     wname(pj)(:ilen(wname(pj))),'=',ww(pj)
 282 c          jj=jj+ilen(wname(pj))+9
 283           write(88,'(bz,2a,e11.5," ",$)')
 284      &     wname(pj)(:ilen(wname(pj))),'=',ww(pj)
 285           jj=jj+ilen(wname(pj))+13
 286         enddo
 287         write (88,'(a,$)') (" ",j=1,79-jj)
 288         write (88,'("&")')
 289       enddo
 290       write (88,'(bz,2(2a,f7.5," "))')
 291      &    "CUTOFF","=",7.0d0,"WCORR4","=",0.0d0
 292       close(88)
 293       if (mod_side .or. mod_side_other) then
 294       if (nparmset.eq.1) then
 295       open(88,
 296      &     file="sc_"//POT(:ilen(POT))//"_opt."//prefix(:ilen(prefix)),
 297      &     status="unknown")
 298       else
 299       open(88,
 300      &     file="sc_"//POT(:ilen(POT))//"_opt."//prefix(:ilen(prefix))//
 301      &     "_par"//liczba,
 302      &     status="unknown")
 303       endif
 304       lprint=.true.
 305       write(88,'(2i5)') ipot,expon
 306       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 307      & ', exponents are ',expon,2*expon
 308       goto (10,20,30,30,40) ipot
 309 C----------------------- LJ potential ---------------------------------
 310    10 write (88,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 311       if (lprint) then
 312         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 313         write (iout,'(a/)') 'The epsilon array:'
 314         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 315         write (iout,'(/a)') 'One-body parameters:'
 316         write (iout,'(a,4x,a)') 'residue','sigma'
 317         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 318       endif
 319       goto 55
 320 C----------------------- LJK potential --------------------------------
 321    20 do i=1,ntyp
 322         write (88,'(4f20.10)')(eps(i,j),j=i,ntyp)
 323         write (iout,*)
 324       enddo
 325       write (88,'(4f20.10)') (sigma0(i),i=1,ntyp)
 326       write (88,*)
 327       write (88,'(4f20.10)') (rr0(i),i=1,ntyp)
 328       write (88,*)
 329       if (lprint) then
 330         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 331         write (iout,'(a/)') 'The epsilon array:'
 332         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 333         write (iout,'(/a)') 'One-body parameters:'
 334         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 335         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 336      &        i=1,ntyp)
 337       endif
 338       goto 55
 339 C---------------------- GB or BP potential -----------------------------
 340    30 do i=1,ntyp
 341         write (88,'(4f20.15)') (eps(i,j),j=i,ntyp)
 342         write (88,*)
 343       enddo
 344       write (88,'(4f20.15)')(sigma0(i),i=1,ntyp)
 345       write (88,*)
 346       write (88,'(4f20.15)')(sigii(i),i=1,ntyp)
 347       write (88,*)
 348       write (88,'(4f20.15)')(chip0(i),i=1,ntyp)
 349       write (88,*)
 350       write (88,'(4f20.15)')(alp(i),i=1,ntyp)
 351       write (88,*)
 352 C For the GB potential convert sigma'**2 into chi'
 353       if (ipot.eq.4) then
 354         do i=1,ntyp
 355           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 356         enddo
 357       endif
 358       if (lprint) then
 359         if (ipot.eq.3) then
 360           write (iout,'(/a/)') 'Parameters of the BP potential:'
 361         else
 362           write (iout,'(/a/)') 'Parameters of the GB potential:'
 363         endif
 364         write (iout,'(a/)') 'The epsilon array:'
 365         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 366         write (iout,'(/a)') 'One-body parameters (asterisk: optimized):'
 367         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 368      &       '    chip  ','    alph  '
 369         do i=1,ntyp
 370         write (iout,'(a3,6x,f10.5,1x,a1,3(f8.5,1x,a1))') restyp(i),
 371      & sigma0(i),aster(mask_sigma(1,i)),sigii(i),aster(mask_sigma(2,i)),
 372      &  chip(i),aster(mask_sigma(3,i)),alp(i),aster(mask_sigma(4,i))
 373         enddo
 374       endif
 375       goto 55
 376 C--------------------- GBV potential -----------------------------------
 377    40 write (88,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 378      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 379      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 380       if (lprint) then
 381         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 382         write (iout,'(a/)') 'The epsilon array:'
 383         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 384         write (iout,'(/a)') 'One-body parameters:'
 385         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 386      &      's||/s_|_^2','    chip  ','    alph  '
 387         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 388      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 389       endif
 390    55 continue
 391 C-----------------------------------------------------------------------
 392       close(88)
 393       endif
 394 #ifdef NEWCORR
 395       if (mod_tor .or. tor_mode.eq.2) then
 396 #else
 397       if (mod_tor) then
 398 #endif
 399         write (iout,'(a)') "Optimized weights of the torsional terms."
 400         do iblock=1,2
 401         do i=-ntortyp+1,ntortyp-1
 402           do j=-ntortyp+1,ntortyp-1
 403             if (mask_tor(0,j,i,iblock).gt.0) then
 404               write (iout,'(i4,2a4,f10.5)') 0,restyp(iloctyp(i)),
 405      &           restyp(iloctyp(j)),weitor(0,j,i,iblock)
 406             endif
 407           enddo
 408         enddo
 409         enddo
 410         if (nparmset.eq.1) then
 411
 412         open(88,file="tor_opt.parm."//prefix(:ilen(prefix)),
 413      &   status="unknown")
 414         write (88,'(i1,7h  ***  ,2a,7h  ***  )') ntortyp,
 415      &   "Optimized torsional parameters ",prefix(:ilen(prefix))
 416         write (88,'(40i2)') (itortyp(i),i=1,ntyp)
 417
 418         if (tor_mode.eq.0) then
 419
 420         write (iout,'(/a/)') 'Torsional constants:'
 421         do iblock=1,2
 422         do i=0,ntortyp-1
 423           do j=0,ntortyp-1
 424             write (iout,*) 'ityp',i,' jtyp',j
 425             write (iout,*) 'Fourier constants'
 426             do k=1,nterm(i,j,iblock)
 427               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 428      &        v2(k,i,j,iblock)
 429             enddo
 430             write (iout,*) 'Lorenz constants'
 431             do k=1,nlor(i,j,iblock)
 432               write (iout,'(3(1pe15.5))')
 433      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 434             enddo
 435           enddo
 436         enddo
 437         enddo
 438
 439         do iblock=1,2
 440         do i=0,ntortyp-1
 441           do j=-ntortyp+1,ntortyp-1
 442             write (88,'(2i2,16h  ************  ,a1,1h-,a1)')
 443      &       nterm(i,j,iblock),nlor(i,j,iblock),onelet(iloctyp(i)),
 444      &       onelet(iloctyp(j))
 445             do k=1,nterm(i,j,iblock)
 446               write (88,'(i6,13x,2(1pe18.5))') k,
 447      &         v1(k,i,j,iblock)*weitor(0,i,j,iblock),v2(k,i,j,
 448      &         iblock)*weitor(0,i,j,iblock)
 449             enddo
 450             do k=1,nlor(i,j,iblock)
 451               write (88,'(i6,13x,3(1pe18.5))')
 452      &         k,vlor1(k,i,j)*weitor(k,i,j,iblock),
 453      &         vlor2(k,i,j)*weitor(k,i,j,iblock),
 454      &         vlor3(k,i,j)*weitor(k,i,j,iblock)
 455             enddo
 456           enddo
 457         enddo
 458         enddo
 459
 460         else
 461 c Print valence-torsional parameters
 462         write (iout,'(a)')
 463      &    "Parameters of the valence-torsional potentials"
 464         do i=-ntortyp+1,ntortyp-1
 465         do j=-ntortyp+1,ntortyp-1
 466         write (iout,'(3a)')"Type ",onelet(iloctyp(i)),onelet(iloctyp(j))
 467         write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
 468         do k=1,nterm_kcc(j,i)
 469           do l=1,nterm_kcc_Tb(j,i)
 470             do ll=1,nterm_kcc_Tb(j,i)
 471                write (iout,'(3i5,2f15.4)')
 472      &            k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
 473             enddo
 474           enddo
 475         enddo
 476         enddo
 477         enddo
 478
 479         do i=-ntortyp+1,ntortyp-1
 480           do j=-ntortyp+1,ntortyp-1
 481             write (88,'(16h  ************  ,a1,1h-,a1)')
 482      &       onelet(iloctyp(i)),onelet(iloctyp(j))
 483             write (88,'(2i5)') nterm_kcc(j,i),nterm_kcc_Tb(j,i)
 484             do k=1,nterm_kcc(j,i)
 485               do l=1,nterm_kcc_Tb(j,i)
 486                 do ll=1,nterm_kcc_Tb(j,i)
 487                   write (88,'(3i5,2(1pe15.5))') k,l-1,ll-1,
 488      &            v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
 489                 enddo
 490               enddo
 491             enddo
 492           enddo
 493         enddo
 494
 495         close(88)
 496
 497         endif
 498
 499         endif
 500       endif
 501       if (mod_sccor) then
 502         write (iout,'(a)') "Optimized weights of the sccor terms."
 503         do i=1,nsccortyp
 504           do j=1,nsccortyp
 505             do l=1,3
 506             if (mask_tor(l,j,i,iblock).gt.0) then
 507               write (iout,'(i4,2a4,f10.5)') l,restyp(isccortyp(j)),
 508      &           restyp(isccortyp(i)),weitor(l,j,i,1)
 509             endif
 510             enddo
 511           enddo
 512         enddo
 513
 514         if (nparmset.eq.1) then
 515
 516         open(88,file="sccor_opt.parm."//prefix(:ilen(prefix)),
 517      &     status="unknown")
 518         write (88,'(i2,7h  ***  ,2a,7h  ***  )') nsccortyp,
 519      & "Optimized local correlation parameters ",prefix(:ilen(prefix))
 520         write (88,'(20i4)') (isccortyp(i),i=1,ntyp)
 521         do l=1,3
 522         do i=1,nsccortyp
 523           do j=1,nsccortyp
 524             write (88,'(2i4)') nterm_sccor(i,j),nlor_sccor(i,j)
 525             write (88,'(12(1h*),1x,a1,1h-,a1)') onelet(i),onelet(j)
 526             do k=1,nterm_sccor(i,j)
 527               write (88,'(i6,13x,2(1pe18.5))') k,
 528      &          v1sccor(k,l,i,j)*weitor(l,i,j,1),
 529      &          v2sccor(k,l,i,j)*weitor(l,i,j,1)
 530             enddo
 531           enddo
 532         enddo
 533         enddo
 534         close(88)
 535
 536         endif
 537
 538       endif
 539 C-----------------------------------------------------------------------
 540 c 7/8/17 AL: Optimization of the bending parameters
 541       if (mod_ang) then
 542         write(iout,*)
 543      &    "Optimized angle potential coefficients (same for D-types)"
 544         do i=0,nthetyp-1
 545           if (mask_ang(i).eq.0) cycle
 546           write (iout,*) "Type: ",onelet(iloctyp(i))
 547           do j=1,nbend_kcc_TB(i)
 548             write (iout,'(i5,f15.5)') j,v1bend_chyb(j,i)
 549           enddo
 550         enddo
 551         open(88,file="theta_opt.parm."//prefix(:ilen(prefix)),
 552      &   status="unknown")
 553         write (88,'(i2)') nthetyp
 554         do i=-nthetyp+1,nthetyp-1
 555           write (88,'(i2,1x,12(1h*),1x,a1)')
 556      &      nbend_kcc_TB(i),onelet(iloctyp(i))
 557           do j=0,nbend_kcc_TB(i)
 558             write (88,'(i5,f20.10)') j,v1bend_chyb(j,i)
 559           enddo
 560         enddo
 561         close(88)
 562       endif
 563 C-----------------------------------------------------------------------
 564       if(mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)then
 565 #ifdef NEWCORR
 566         write (iout,'(2a)')
 567      &   "Fourier coefficients of new (angle-dependent) cumulants",
 568      &   " (asterisk: optimized)"
 569         do i=0,nloctyp-1
 570           write (iout,*) "Type: ",onelet(iloctyp(i))
 571           write (iout,*) "Coefficients of the expansion of B1"
 572           do j=1,2
 573             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
 574           enddo
 575           write (iout,*) "Coefficients of the expansion of B2"
 576           do j=1,2
 577             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
 578           enddo
 579           write (iout,*) "Coefficients of the expansion of C"
 580           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
 581           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
 582           write (iout,*) "Coefficients of the expansion of D"
 583           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
 584           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
 585           write (iout,*) "Coefficients of the expansion of E"
 586           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
 587           do j=1,2
 588             do k=1,2
 589               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
 590             enddo
 591           enddo
 592         enddo
 593
 594         IF (SPLIT_FOURIERTOR) THEN
 595
 596         write (iout,'(2a)')
 597      &   "Fourier coefficients of new (angle-dependent) cumulants",
 598      &   " for torsional potentials (asterisk: optimized)"
 599         do i=0,nloctyp-1
 600           write (iout,*) "Type: ",onelet(iloctyp(i))
 601           write (iout,*) "Coefficients of the expansion of B1"
 602           do j=1,2
 603             write (iout,'(3hB1(,i1,1h),3f10.5)')
 604      &         j,(bnew1tor(k,j,i),k=1,3)
 605           enddo
 606           write (iout,*) "Coefficients of the expansion of B2"
 607           do j=1,2
 608             write (iout,'(3hB2(,i1,1h),3f10.5)')
 609      &         j,(bnew2tor(k,j,i),k=1,3)
 610           enddo
 611           write (iout,*) "Coefficients of the expansion of C"
 612           write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
 613           write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
 614           write (iout,*) "Coefficients of the expansion of D"
 615           write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
 616           write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
 617           write (iout,*) "Coefficients of the expansion of E"
 618           write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
 619           do j=1,2
 620             do k=1,2
 621               write (iout,'(1hE,2i1,2f10.5)')
 622      &         j,k,(eenewtor(l,j,k,i),l=1,2)
 623             enddo
 624           enddo
 625         enddo
 626
 627         ENDIF
 628
 629 #else
 630         write (iout,'(2a)')
 631      &   "Fourier coefficients of old angle-independent cumulants",
 632      &   " (asterisk: optimized)"
 633         do i=0,nloctyp-1
 634           write (iout,*) 'Type ',restyp(iloctyp(i))
 635             write (iout,'(a,i2,a,f10.5,1x,a)')
 636      &    ('b(',j,')=',b(j,i),aster(mysign(mask_fourier(j,i))),j=1,13)
 637         enddo
 638 #endif
 639 c Write a file with parameters for UNRES
 640         if (nparmset.eq.1) then
 641         open(88,file="fourier_opt.parm."//prefix(:ilen(prefix)),
 642      &     status="unknown")
 643         else
 644         open(88,file="fourier_opt.parm."//prefix(:ilen(prefix))//
 645      &     "_par"//liczba,
 646      &     status="unknown")
 647         endif
 648 #ifdef NEWCORR
 649         IF (SPLIT_FOURIERTOR) THEN
 650           write (88,'(i5,5x,a)') -nloctyp,
 651      &                         "# Number of local interaction types"
 652         ELSE
 653           write (88,'(i5,5x,a)') nloctyp,
 654      &                         "# Number of local interaction types"
 655         ENDIF
 656         write (88,'(40i2)') (itype2loc(i),i=1,ntyp)
 657         write (88,'(20i4)') (iloctyp(i),i=0,nloctyp)
 658         do i=0,nloctyp-1
 659           write (88,'(a)') restyp(iloctyp(i))
 660           do ii=1,3
 661             do j=1,2
 662               write (88,'(f20.10,4h  b1,i1,1h(,i1,1h))')
 663      &         bnew1(ii,j,i),j,ii
 664             enddo
 665           enddo
 666           do ii=1,3
 667             do j=1,2
 668               write (88,'(f20.10,4h  b2,i1,1h(,i1,1h))')
 669      &         bnew2(ii,j,i),j,ii
 670             enddo
 671           enddo
 672           do ii=1,3
 673             do j=1,2
 674               write (88,'(f20.10,4h  c1,i1,1h(,i1,1h))')
 675      &         2*ccnew(ii,j,i),j,ii
 676             enddo
 677           enddo
 678           do ii=1,3
 679             do j=1,2
 680               write (88,'(f20.10,4h  d1,i1,1h(,i1,1h))')
 681      &         2*ddnew(ii,j,i),j,ii
 682             enddo
 683           enddo
 684           do ii=1,2
 685             do j=1,2
 686               do k=1,2
 687               write (88,'(f20.10,3h  e,2i1,1h(,i1,1h))')
 688      &         eenew(ii,j,k,i),j,k,ii
 689               enddo
 690             enddo
 691           enddo
 692           do ii=1,3
 693               write (88,'(f20.10,5h  e0(,i1,1h))')
 694      &         e0new(ii,i),ii
 695           enddo
 696         enddo
 697
 698         IF (SPLIT_FOURIERTOR) THEN
 699
 700         do i=0,nloctyp-1
 701           write (88,'(a)') restyp(iloctyp(i))
 702           do ii=1,3
 703             do j=1,2
 704               write (88,'(f20.10,4h  b1,i1,1h(,i1,1h))')
 705      &         bnew1tor(ii,j,i),j,ii
 706             enddo
 707           enddo
 708           do ii=1,3
 709             do j=1,2
 710               write (88,'(f20.10,4h  b2,i1,1h(,i1,1h))')
 711      &         bnew2tor(ii,j,i),j,ii
 712             enddo
 713           enddo
 714           do ii=1,3
 715             do j=1,2
 716               write (88,'(f20.10,4h  c1,i1,1h(,i1,1h))')
 717      &         2*ccnewtor(ii,j,i),j,ii
 718             enddo
 719           enddo
 720           do ii=1,3
 721             do j=1,2
 722               write (88,'(f20.10,4h  d1,i1,1h(,i1,1h))')
 723      &         2*ddnewtor(ii,j,i),j,ii
 724             enddo
 725           enddo
 726           do ii=1,2
 727             do j=1,2
 728               do k=1,2
 729               write (88,'(f20.10,3h  e,2i1,1h(,i1,1h))')
 730      &         eenewtor(ii,j,k,i),j,k,ii
 731               enddo
 732             enddo
 733           enddo
 734           do ii=1,3
 735               write (88,'(f20.10,5h  e0(,i1,1h))')
 736      &         e0newtor(ii,i),ii
 737           enddo
 738         enddo
 739
 740         ENDIF
 741
 742 #else
 743         write (88,'(i5,5x,a)') nloctyp,
 744      &                         "# Number of local interaction types"
 745         do i=1,nloctyp
 746           write (88,'(a)') restyp(iloctyp(i))
 747           do j=1,13
 748             write (88,'(f20.15)') b(j,i)
 749           enddo
 750         enddo
 751 #endif
 752         close(88)
 753       endif
 754       if (mod_elec) then
 755                 write (iout,'(/a)')
 756      &  'Electrostatic interaction constants (asterisk: optimized):'
 757         write (iout,'(1x,a,1x,a,7x,a,15x,a,15x,a,13x,a)')
 758      &            'IT','JT','EPP','RPP','ELPP6','ELPP3'
 759         do i=1,2
 760           do j=1,2
 761             write(iout,'(2i3,4(1pe15.4,1x,a1,1x))')i,j,
 762      &       epp(i,j),aster(mask_elec(i,j,1)),
 763      &       rpp(i,j),aster(mask_elec(i,j,2)),
 764      &       elpp6(i,j),aster(mask_elec(i,j,3)),
 765      &       elpp3(i,j),aster(mask_elec(i,j,4))
 766           enddo
 767         enddo
 768         write (iout,*)
 769         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 770      &            'IT','JT','APP','BPP','AEL6','AEL3'
 771         do i=1,2
 772           do j=1,2
 773             write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 774      &                    ael6(i,j),ael3(i,j)
 775           enddo
 776         enddo
 777 C Write electrostatic interaction parameters to a file
 778         if (nparmset.eq.1) then
 779         open(88,file="electr_opt.parm."//prefix(:ilen(prefix)),
 780      &   status="unknown")
 781         else
 782         open(88,file="electr_opt.parm."//prefix(:ilen(prefix))//
 783      &   "_par"//liczba,
 784      &   status="unknown")
 785         endif
 786         write(88,'(4f15.10,5x,a)')((epp(i,j),j=1,2),i=1,2),"! EPP"
 787         write(88,'(4f15.10,5x,a)')((rpp(i,j),j=1,2),i=1,2),"! RPP"
 788         write(88,'(4f15.10,5x,a)')((elpp6(i,j),j=1,2),i=1,2),"! ELPP6"
 789         write(88,'(4f15.10,5x,a)')((elpp3(i,j),j=1,2),i=1,2),"! ELPP3"
 790         close(88)
 791       endif
 792       if (mod_scp) then
 793         write (iout,*)
 794         write (iout,*)
 795      &  "Parameters of SC-p interactions (star: optimized):"
 796         write (iout,'(t14,a,t34,a)') "TYPE 1","TYPE 2"
 797         write (iout,'(8a)') "SC TYPE"," EPS_SCP  ","   R_SCP  ",
 798      &    " EPS_SCP  ","   R_SCP  "
 799         do i=1,20
 800           ast11 = aster(mask_scp(i,1,1))
 801           ast12 = aster(mask_scp(i,1,2))
 802           ast21 = aster(mask_scp(i,2,1))
 803           ast22 = aster(mask_scp(i,2,2))
 804           if (mask_scp(i,1,1).eq.0) ast11 = aster(mask_scp(0,1,1))
 805           if (mask_scp(i,1,2).eq.0) ast12 = aster(mask_scp(0,1,2))
 806           if (mask_scp(i,2,1).eq.0) ast21 = aster(mask_scp(0,2,1))
 807           if (mask_scp(i,2,2).eq.0) ast22 = aster(mask_scp(0,2,2))
 808           write (iout,'(2x,A3,2x,4(f8.3,1x,a1))') restyp(i),
 809      &     eps_scp(i,1),ast11,rscp(i,1),ast12,eps_scp(i,2),ast21,
 810      &     rscp(i,2),ast22
 811         enddo
 812 C Write SCp parameters to a file
 813         if (nparmset.eq.1) then
 814         open(88,file="scp_opt.parm."//prefix(:ilen(prefix)),
 815      &     status="unknown")
 816         else
 817         open(88,file="scp_opt.parm."//prefix(:ilen(prefix))//
 818      &     "_par"//liczba,
 819      &     status="unknown")
 820         endif
 821         do i=1,ntyp
 822           write(88,'(4f15.10,5x,a)') eps_scp(i,1),rscp(i,1),
 823      &      eps_scp(i,2),rscp(i,2),restyp(i)
 824         enddo
 825         close(88)
 826       endif
 827 50    format(bz,15(f8.5,1x,a1))
 828 60    format(bz,100f10.5)
 829       end
 830 c------------------------------------------------------------------------------
 831       subroutine write_zscore
 832       implicit none
 833       include "DIMENSIONS"
 834       include "DIMENSIONS.ZSCOPT"
 835       include "COMMON.ENERGIES"
 836       include "COMMON.CLASSES"
 837       include "COMMON.PROTNAME"
 838       include "COMMON.VMCPAR"
 839       include "COMMON.IOUNITS"
 840       include "COMMON.WEIGHTS"
 841       include "COMMON.WEIGHTDER"
 842       include "COMMON.COMPAR"
 843       include "COMMON.OPTIM"
 844       include "COMMON.THERMAL"
 845       integer i,j,k,iprot,ib
 846       character*6 namnat(5) /"angrms","Q","rmsd","rgy","sign"/
 847       integer ilen
 848       external ilen
 849       do iprot=1,nprot
 850
 851         write (iout,*)
 852         write (iout,'(a,2x,a)') "Protein",
 853      &    protname(iprot)(:ilen(protname(iprot)))
 854
 855         write (iout,*)
 856         write (iout,'(a)')"Maximum likelihood."
 857         write (iout,*)
 858         do ib=1,nbeta(1,iprot)
 859           write (iout,'(f8.1,f10.5)')
 860      &      1.0d0/(Rgas*betaT(ib,1,iprot)),sumlik(ib,iprot)
 861         enddo
 862         write (iout,*)
 863
 864         write (iout,*)
 865         write (iout,'(a)')"Specific heat."
 866         write (iout,*)
 867         write (iout,'(a,f10.5)') "Baseline",heatbase(iprot)
 868         do ib=1,nbeta(2,iprot)
 869           write (iout,'(f8.1,3f10.5)')
 870      &      1.0d0/(Rgas*betaT(ib,2,iprot)),
 871      &      zvtot(ib,iprot),target_cv(ib,iprot),weiCv(ib,iprot)
 872         enddo
 873
 874         if (natlike(iprot).gt.0) write (iout,'(a)')"Nativelikeness."
 875         do i=1,natlike(iprot)
 876           write (iout,*) "Property",i,numnat(i,iprot)
 877           if (natdim(i,iprot).eq.1)  then
 878             do ib=1,nbeta(i+2,iprot)
 879               write (iout,'(f7.2,2f10.5)')
 880      &         1.0d0/(Rgas*betaT(ib,i+2,iprot)),nuave(1,ib,i,iprot),
 881      &         nuexp(1,ib,i,iprot),weinat(1,ib,i,iprot)
 882             enddo
 883           else
 884             do ib=1,nbeta(i+2,iprot)
 885               write (iout,'("Temperature",f7.2," NATDIM",i5)')
 886      &        1.0d0/(Rgas*betaT(ib,i+2,iprot)),natdim(i,iprot)
 887               do k=1,natdim(i,iprot)
 888               write (iout,'(3f10.5)') nuave(k,ib,i,iprot),
 889      &         nuexp(k,ib,i,iprot),weinat(k,ib,i,iprot)
 890               enddo
 891             enddo
 892           endif
 893         enddo
 894
 895       enddo
 896       write (iout,*)
 897       if (torsion_pmf.or.turn_pmf.or.eello_pmf)
 898      &  write (iout,*) "PMF contribution to the target function",fpmf
 899       if (pdbpmf) write (iout,*)
 900      & "PDB angle-torsional statistics contribution",fpdb
 901       write (iout,*)
 902       return
 903       end
 904 c------------------------------------------------------------------------------
 905       subroutine write_prot(ii,przed)
 906       implicit none
 907       include "DIMENSIONS.ZSCOPT"
 908       include "COMMON.ENERGIES"
 909       include "COMMON.CLASSES"
 910       include "COMMON.PROTNAME"
 911       include "COMMON.VMCPAR"
 912       include "COMMON.IOUNITS"
 913       include "COMMON.COMPAR"
 914       integer i,ii,j
 915       integer ilen
 916       external ilen
 917       character*(*) przed
 918       character*64 nazwa
 919       nazwa=przed(:ilen(przed))//'.'//prefix(:ilen(prefix))
 920      &   //'.'//protname(ii)
 921       open(unit=istat,file=nazwa)
 922       do i=1,ntot_work(ii)
 923         write(istat,'(i10,3f15.3,20e15.5)')
 924      &    i,e_total(i,ii),eini(i,ii),entfac(i,ii),
 925      &    (Ptab(i,j,ii),j=1,nbeta(1,ii))
 926       enddo
 927       close(istat)
 928       return
 929       end
 930 c------------------------------------------------------------------------------
 931       integer function mysign(i)
 932       integer  i
 933       if (i.gt.0) then
 934         mysign=1
 935       else
 936         mysign=0
 937       endif
 938       return
 939       end