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