source/wham/src-old/enecalc1.F

   1       subroutine enecalc(islice,*)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5       include "DIMENSIONS.FREE"
   6 #ifdef MPI
   7       include "mpif.h"
   8       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
   9       include "COMMON.MPI"
  10 #endif
  11       include "COMMON.CHAIN"
  12       include "COMMON.IOUNITS"
  13       include "COMMON.PROTFILES"
  14       include "COMMON.NAMES"
  15       include "COMMON.VAR"
  16       include "COMMON.SBRIDGE"
  17       include "COMMON.GEO"
  18       include "COMMON.FFIELD"
  19       include "COMMON.ENEPS"
  20       include "COMMON.LOCAL"
  21       include "COMMON.WEIGHTS"
  22       include "COMMON.INTERACT"
  23       include "COMMON.FREE"
  24       include "COMMON.ENERGIES"
  25       include "COMMON.CONTROL"
  26       include "COMMON.TORCNSTR"
  27       character*64 nazwa
  28       character*80 bxname
  29       character*3 liczba
  30       double precision qwolynes
  31       external qwolynes
  32       integer errmsg_count,maxerrmsg_count /100/
  33       double precision rmsnat,gyrate
  34       external rmsnat,gyrate
  35       double precision tole /1.0d-1/
  36       integer i,itj,ii,iii,j,k,l,licz,scme,itmp
  37       integer ires
  38       integer ir,ib,ipar,iparm
  39       integer iscor,islice
  40       real*4 csingle(3,maxres2)
  41       double precision energ
  42       integer ilen,iroof
  43       external ilen,iroof
  44       double precision energia(0:max_ene),rmsdev,efree,eini
  45       double precision fT(6),quot,quotl,kfacl,kfac /2.4d0/,T0 /3.0d2/
  46       double precision tt
  47       integer snk_p(MaxR,MaxT_h,Max_parm)
  48       logical lerr
  49       character*64 bprotfile_temp
  50       integer scount_t(0:maxprocs-1)
  51       call opentmp(islice,ientout,bprotfile_temp)
  52       iii=0
  53       ii=0
  54       errmsg_count=0
  55       write (iout,*) "enecalc: nparmset ",nparmset
  56 #ifdef MPI
  57       do iparm=1,nParmSet
  58         do ib=1,nT_h(iparm)
  59           do i=1,nR(ib,iparm)
  60             snk_p(i,ib,iparm)=0
  61           enddo
  62         enddo
  63       enddo
  64       do i=indstart(me1),indend(me1)
  65 #else
  66       do iparm=1,nParmSet
  67         do ib=1,nT_h(iparm)
  68           do i=1,nR(ib,iparm)
  69             snk(i,ib,iparm)=0
  70           enddo
  71         enddo
  72       enddo
  73       do i=1,ntot
  74 #endif
  75         read(ientout,rec=i,err=101)
  76      &    ((csingle(l,k),l=1,3),k=1,nres),
  77      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
  78      &    nss,(ihpb(k),jhpb(k),k=1,nss),
  79      &    eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
  80 #ifdef DEBUG
  81         write (iout,*) "homol_nset",homol_nset,
  82      &    "  i",i," iR",iR," ib",ib," iset",iset
  83 #endif
  84         if (homol_nset.gt.1) iset=iR
  85 #ifdef DEBUG
  86         write (iout,*) "homol_nset",homol_nset,
  87      &    "  i",i," iR",iR," ib",ib," iset",iset
  88 #endif
  89 cc       write(iout,*), 'NAWEJ',i,eini
  90          if (indpdb.gt.0) then
  91            do k=1,nres
  92              do l=1,3
  93                c(l,k)=csingle(l,k)
  94              enddo
  95            enddo
  96            do k=nnt,nct
  97              do l=1,3
  98                c(l,k+nres)=csingle(l,k+nres)
  99              enddo
 100            enddo
 101            q(nQ+1,iii+1)=rmsnat(iii+1)
 102          endif
 103          q(nQ+2,iii+1)=gyrate(iii+1)
 104 c        fT=T0*beta_h(ib,ipar)*1.987D-3
 105 c        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
 106         if (rescale_mode.eq.1) then
 107           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 108 #if defined(FUNCTH)
 109           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 110           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 111 #elif defined(FUNCT)
 112           ft(6)=quot
 113 #else
 114           ft(6)=1.0d0
 115 #endif
 116           quotl=1.0d0
 117           kfacl=1.0d0
 118           do l=1,5
 119             quotl=quotl*quot
 120             kfacl=kfacl*kfac
 121             fT(l)=kfacl/(kfacl-1.0d0+quotl)
 122           enddo
 123         else if (rescale_mode.eq.2) then
 124           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 125 #if defined(FUNCTH)
 126           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 127           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 128 #elif defined(FUNCT)
 129           ft(6)=quot
 130 #else
 131           ft(6)=1.0d0
 132 #endif
 133           quotl=1.0d0
 134           do l=1,5
 135             quotl=quotl*quot
 136             fT(l)=1.12692801104297249644d0/
 137      &         dlog(dexp(quotl)+dexp(-quotl))
 138           enddo
 139         else if (rescale_mode.eq.0) then
 140           do l=1,5
 141             fT(l)=1.0d0
 142           enddo
 143         else
 144           write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",
 145      &     rescale_mode
 146           call flush(iout)
 147           return1
 148         endif
 149
 150 c        write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
 151 c     &   " kfac",kfac,"quot",quot," fT",fT
 152         do j=1,2*nres
 153           do k=1,3
 154             c(k,j)=csingle(k,j)
 155           enddo
 156         enddo
 157         call int_from_cart1(.false.)
 158         ii=ii+1
 159         do iparm=1,nparmset
 160
 161         call restore_parm(iparm)
 162 #ifdef DEBUG
 163             write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 164      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 165      &      wtor_d,wsccor,wbond
 166 #endif
 167         call etotal(energia(0),fT)
 168         if (constr_homology) energia(0)=energia(0)+
 169      &   waga_homology(iset)*energia(22)
 170         write (iout,*) "constr_homology",constr_homology," iset",iset,
 171      &   " waga_homology",waga_homology(iset)
 172 #ifdef DEBUG
 173         write (iout,*) "Conformation",i
 174         call enerprint(energia(0),fT)
 175 c        write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21)
 176 c        write (iout,*) "ftors",ftors
 177 c      write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 178 c     & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 179 c      do ires=1,nres
 180 c        write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 181 c     &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 182 c     &    (c(j,ires+nres),j=1,3)
 183 c      enddo
 184         call intout
 185 #endif
 186         if (energia(0).ge.1.0d20) then
 187           write (iout,*) "NaNs detected in some of the energy",
 188      &     " components for conformation",ii+1
 189           write (iout,*) "The Cartesian geometry is:"
 190           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 191           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 192           write (iout,*) "The internal geometry is:"
 193 c          call intout
 194 c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 195           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 196           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 197           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 198           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 199           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 200           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 201           write (iout,*) "The components of the energy are:"
 202           call enerprint(energia(0),fT)
 203           write (iout,*)
 204      &      "This conformation WILL NOT be added to the database."
 205           call flush(iout)
 206           goto 121
 207         else
 208 #ifdef DEBUG
 209           if (ipar.eq.iparm) write (iout,*) i,iparm,
 210      &      1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
 211 #endif
 212           if (ipar.eq.iparm .and. einicheck.gt.0 .and.
 213      &      dabs(eini-energia(0)).gt.tole) then
 214             if (errmsg_count.le.maxerrmsg_count) then
 215               write (iout,'(2a,2e15.5,a,2i8,a,f8.1)')
 216      &         "Warning: energy differs remarkably from ",
 217      &         " the value read in: ",energia(0),eini," point",
 218      &         iii+1,indstart(me1)+iii," T",
 219      &         1.0d0/(1.987D-3*beta_h(ib,ipar))
 220 #ifdef DEBUG
 221              call enerprint(energia(0),fT)
 222              write (iout,'(4f10.5,2i5)') 0.0,energia(0),0.0,
 223      &       1.0d0/(beta_h(ib,ipar)*1.987D-3),
 224      &       0,0
 225              write(iout,'(8f10.5)')
 226      &       ((c(l,k),l=1,3),k=1,nres),
 227      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 228              itmp=ipdb
 229              ipdb=iout
 230              call pdbout(iii+1,beta_h(ib,ipar),
 231      &                   eini,energia(0),0.0d0,rmsdev)
 232              write (iout,*)
 233              ipdb=itmp
 234 #endif
 235               errmsg_count=errmsg_count+1
 236               if (errmsg_count.gt.maxerrmsg_count)
 237      &          write (iout,*) "Too many warning messages"
 238               if (einicheck.gt.1) then
 239                 write (iout,*) "Calculation stopped."
 240                 call flush(iout)
 241 #ifdef MPI
 242                 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
 243 #endif
 244                 call flush(iout)
 245                 return1
 246               endif
 247             endif
 248           endif
 249           potE(iii+1,iparm)=energia(0)
 250           do k=1,max_ene
 251             enetb(k,iii+1,iparm)=energia(k)
 252           enddo
 253 c          write(iout,*) "iCHUJ TU STRZELI",i,iii,entfac(i)
 254 c          call enerprint(energia(0),fT)
 255 #ifdef DEBUG
 256           write (iout,'(2i5,f10.1,3e15.5)') i,iii,
 257      &     1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
 258           call enerprint(energia(0),fT)
 259           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 260           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 261           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 262           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 263           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 264           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 265           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 266           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 267           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 268           write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
 269           write (iout,'(f10.5,i10)') rmsdev,iscor
 270           call enerprint(energia(0),fT)
 271         write(liczba,'(bz,i3.3)') me
 272         nazwa="test"//liczba//".pdb"
 273         write (iout,*) "pdb file",nazwa
 274         open (ipdb,file=nazwa,position="append")
 275         call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 276         close(ipdb)
 277 #endif
 278         endif
 279
 280         enddo ! iparm
 281
 282         iii=iii+1
 283         if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
 284         write (ientout,rec=iii)
 285      &   ((csingle(l,k),l=1,3),k=1,nres),
 286      &   ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 287      &   nss,(ihpb(k),jhpb(k),k=1,nss),
 288      &   potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
 289 c        write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
 290 #ifdef MPI
 291         if (separate_parset) then
 292           snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
 293         else
 294           snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
 295         endif
 296 c        write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
 297 c     &   " snk",snk_p(iR,ib,ipar)
 298 #else
 299         snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
 300 #endif
 301   121   continue
 302       enddo
 303 #ifdef MPI
 304       scount(me)=iii
 305       write (iout,*) "Me",me," scount",scount(me)
 306       call flush(iout)
 307 c  Master gathers updated numbers of conformations written by all procs.
 308       scme = scount(me)
 309       call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_t(0), 1,
 310      &  MPI_INTEGER, WHAM_COMM, IERROR)
 311       do k=0,nprocs-1
 312         scount(k) = scount_t(k)
 313       enddo
 314       indstart(0)=1
 315       indend(0)=scount(0)
 316       do i=1, Nprocs-1
 317         indstart(i)=indend(i-1)+1
 318         indend(i)=indstart(i)+scount(i)-1
 319       enddo
 320       write (iout,*)
 321       write (iout,*) "Revised conformation counts"
 322       do i=0,nprocs1-1
 323         write (iout,'(a,i5,a,i7,a,i7,a,i7)')
 324      &    "Processor",i," indstart",indstart(i),
 325      &    " indend",indend(i)," count",scount(i)
 326       enddo
 327       call flush(iout)
 328       call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),
 329      &  MaxR*MaxT_h*nParmSet,
 330      &  MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
 331 #endif
 332       stot(islice)=0
 333       do iparm=1,nParmSet
 334         do ib=1,nT_h(iparm)
 335           do i=1,nR(ib,iparm)
 336             stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
 337           enddo
 338         enddo
 339       enddo
 340       write (iout,*) "Revised SNK"
 341       do iparm=1,nParmSet
 342         do ib=1,nT_h(iparm)
 343           write (iout,'("Param",i3," Temp",f6.1,3x,32i8)')
 344      &     iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),
 345      &     (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
 346           write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
 347         enddo
 348       enddo
 349       write (iout,'("Total",i10)') stot(islice)
 350       call flush(iout)
 351       return
 352   101 write (iout,*) "Error in scratchfile."
 353       call flush(iout)
 354       return1
 355       end
 356 c------------------------------------------------------------------------------
 357       subroutine write_dbase(islice,*)
 358       implicit none
 359       include "DIMENSIONS"
 360       include "DIMENSIONS.ZSCOPT"
 361       include "DIMENSIONS.FREE"
 362       include "DIMENSIONS.COMPAR"
 363 #ifdef MPI
 364       include "mpif.h"
 365       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 366       include "COMMON.MPI"
 367 #endif
 368       include "COMMON.CONTROL"
 369       include "COMMON.CHAIN"
 370       include "COMMON.IOUNITS"
 371       include "COMMON.PROTFILES"
 372       include "COMMON.NAMES"
 373       include "COMMON.VAR"
 374       include "COMMON.SBRIDGE"
 375       include "COMMON.GEO"
 376       include "COMMON.FFIELD"
 377       include "COMMON.ENEPS"
 378       include "COMMON.LOCAL"
 379       include "COMMON.WEIGHTS"
 380       include "COMMON.INTERACT"
 381       include "COMMON.FREE"
 382       include "COMMON.ENERGIES"
 383       include "COMMON.COMPAR"
 384       include "COMMON.PROT"
 385       character*64 nazwa
 386       character*80 bxname,cxname
 387       character*64 bprotfile_temp
 388       character*3 liczba,licz
 389       character*2 licz2
 390       integer i,itj,ii,iii,j,k,l
 391       integer ixdrf,iret
 392       integer iscor,islice
 393       double precision rmsdev,efree,eini
 394       real*4 csingle(3,maxres2)
 395       double precision energ
 396       integer ilen,iroof
 397       external ilen,iroof
 398       integer ir,ib,iparm
 399       write (licz2,'(bz,i2.2)') islice
 400       call opentmp(islice,ientout,bprotfile_temp)
 401       write (iout,*) "bprotfile_temp ",bprotfile_temp
 402       call flush(iout)
 403       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0
 404      &   .and. ensembles.eq.0) then
 405         close(ientout,status="delete")
 406         return
 407       endif
 408 #ifdef MPI
 409       write (liczba,'(bz,i3.3)') me
 410       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 411         if (.not.separate_parset) then
 412           bxname = prefix(:ilen(prefix))//liczba//".bx"
 413         else
 414           write (licz,'(bz,i3.3)') myparm
 415           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 416         endif
 417         open (ientin,file=bxname,status="unknown",
 418      &    form="unformatted",access="direct",recl=lenrec1)
 419       endif
 420 #else
 421       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 422         if (nslice.eq.1) then
 423           bxname = prefix(:ilen(prefix))//".bx"
 424         else
 425           bxname = prefix(:ilen(prefix))//
 426      &           "_slice_"//licz2//".bx"
 427         endif
 428         open (ientin,file=bxname,status="unknown",
 429      &    form="unformatted",access="direct",recl=lenrec1)
 430         write (iout,*) "Calculating energies; writing geometry",
 431      & " and energy components to ",bxname(:ilen(bxname))
 432       endif
 433 #if (defined(AIX) && !defined(JUBL))
 434         call xdrfopen_(ixdrf,cxname, "w", iret)
 435 #else
 436         call xdrfopen(ixdrf,cxname, "w", iret)
 437 #endif
 438         if (iret.eq.0) then
 439           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 440           cxfile=.fale.
 441         endif
 442       endif
 443 #endif
 444       if (indpdb.gt.0) then
 445         if (nslice.eq.1) then
 446 #ifdef MPI
 447          if (.not.separate_parset) then
 448            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 449      &       //liczba//'.stat'
 450          else
 451            write (licz,'(bz,i3.3)') myparm
 452            statname=prefix(:ilen(prefix))//'_par'//licz//'_'//
 453      &      pot(:ilen(pot))//liczba//'.stat'
 454          endif
 455
 456 #else
 457           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
 458 #endif
 459         else
 460 #ifdef MPI
 461          if (.not.separate_parset) then
 462           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 463      &      "_slice_"//licz2//liczba//'.stat'
 464          else
 465           write (licz,'(bz,i3.3)') myparm
 466           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 467      &      '_par'//licz//"_slice_"//licz2//liczba//'.stat'
 468          endif
 469 #else
 470           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 471      &      //"_slice_"//licz2//'.stat'
 472 #endif
 473         endif
 474         open(istat,file=statname,status="unknown")
 475       endif
 476
 477 #ifdef MPI
 478       do i=1,scount(me)
 479 #else
 480       do i=1,ntot(islice)
 481 #endif
 482 cc        if (dyn_ss) then
 483 cc        read(ientout,rec=i,err=101)
 484 cc     &    ((csingle(l,k),l=1,3),k=1,nres),
 485 cc     &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 486 cc     &    nss,(idssb(k),jdssb(k),k=1,nss),
 487 cc     &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
 488 cc        idssb(k)=idssb(k)-nres
 489 cc        jdssb(k)=jdssb(k)-nres
 490 cc        else
 491         read(ientout,rec=i,err=101)
 492      &    ((csingle(l,k),l=1,3),k=1,nres),
 493      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 494      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 495      &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
 496 cc         endif
 497 cc        write (iout,*) 'CC', iR,ib,iparm,eini,efree
 498         do j=1,2*nres
 499           do k=1,3
 500             c(k,j)=csingle(k,j)
 501           enddo
 502         enddo
 503         call int_from_cart1(.false.)
 504         iscore=0
 505         if (indpdb.gt.0) then
 506           call conf_compar(i,.false.,.true.)
 507         endif
 508 c        if (dyn_ss) then
 509         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i)
 510      &    ((csingle(l,k),l=1,3),k=1,nres),
 511      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 512      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 513 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
 514      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
 515 c        else
 516           if (bxfile .or.cxfile .or. ensembles.gt.0) write
 517      &     (ientin,rec=i)
 518      &    ((csingle(l,k),l=1,3),k=1,nres),
 519      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 520      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 521 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
 522      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
 523 c        endif
 524 c        write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
 525 #ifndef MPI
 526         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),
 527      &    -entfac(i),rms_nat,iscore)
 528 #endif
 529       enddo
 530       close(ientout,status="delete")
 531       close(istat)
 532       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
 533 #ifdef MPI
 534       call MPI_Barrier(WHAM_COMM,IERROR)
 535       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile
 536      &   .and. ensembles.eq.0) return
 537       write (iout,*)
 538       if (bxfile .or. ensembles.gt.0) then
 539         if (nslice.eq.1) then
 540           if (.not.separate_parset) then
 541             bxname = prefix(:ilen(prefix))//".bx"
 542           else
 543             write (licz,'(bz,i3.3)') myparm
 544             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
 545           endif
 546         else
 547           if (.not.separate_parset) then
 548             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
 549           else
 550             write (licz,'(bz,i3.3)') myparm
 551             bxname = prefix(:ilen(prefix))//"par_"//licz//
 552      &        "_slice_"//licz2//".bx"
 553           endif
 554         endif
 555         open (ientout,file=bxname,status="unknown",
 556      &      form="unformatted",access="direct",recl=lenrec1)
 557         write (iout,*) "Master is creating binary database ",
 558      &   bxname(:ilen(bxname))
 559       endif
 560       if (cxfile) then
 561         if (nslice.eq.1) then
 562           if (.not.separate_parset) then
 563             cxname = prefix(:ilen(prefix))//".cx"
 564           else
 565             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
 566           endif
 567         else
 568           if (.not.separate_parset) then
 569             cxname = prefix(:ilen(prefix))//
 570      &             "_slice_"//licz2//".cx"
 571           else
 572             cxname = prefix(:ilen(prefix))//"_par"//licz//
 573      &             "_slice_"//licz2//".cx"
 574           endif
 575         endif
 576 #if (defined(AIX) && !defined(JUBL))
 577         call xdrfopen_(ixdrf,cxname, "w", iret)
 578 #else
 579         call xdrfopen(ixdrf,cxname, "w", iret)
 580 #endif
 581         if (iret.eq.0) then
 582           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 583           cxfile=.false.
 584         endif
 585       endif
 586       do j=0,nprocs-1
 587         write (liczba,'(bz,i3.3)') j
 588         if (separate_parset) then
 589           write (licz,'(bz,i3.3)') myparm
 590           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 591         else
 592           bxname = prefix(:ilen(prefix))//liczba//".bx"
 593         endif
 594         open (ientin,file=bxname,status="unknown",
 595      &    form="unformatted",access="direct",recl=lenrec1)
 596         write (iout,*) "Master is reading conformations from ",
 597      &   bxname(:ilen(bxname))
 598         iii = 0
 599 c        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
 600 c        call flush(iout)
 601         do i=indstart(j),indend(j)
 602           iii = iii+1
 603 cc          if (dyn_ss) then
 604 cc          read(ientin,rec=iii,err=101)
 605 cc     &      ((csingle(l,k),l=1,3),k=1,nres),
 606 cc     &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 607 cc     &      nss,(idssb(k),jdssb(k),k=1,nss),
 608 cc    &      eini,efree,rmsdev,iscor
 609 cc          idssb(k)=idssb(k)-nres
 610 cc          jdssb(k)=jdssb(k)-nres
 611 cc          else
 612           read(ientin,rec=iii,err=101)
 613      &      ((csingle(l,k),l=1,3),k=1,nres),
 614      &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 615      &      nss,(ihpb(k),jhpb(k),k=1,nss),
 616      &      eini,efree,rmsdev,iscor
 617 cc          endif
 618           if (bxfile .or. ensembles.gt.0) then
 619 cc          if (dyn_ss) then
 620 cc            write (ientout,rec=i)
 621 cc     &        ((csingle(l,k),l=1,3),k=1,nres),
 622 cc     &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 623 cc     &        nss,(idssb(k)+nres,jdssb(k)+nres,k=1,nss),
 624 cc    &        eini,efree,rmsdev,iscor
 625 cc           else
 626                         write (ientout,rec=i)
 627      &        ((csingle(l,k),l=1,3),k=1,nres),
 628      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 629      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 630      &        eini,efree,rmsdev,iscor
 631 cc           write(iout,*) "W poszukiwaniu zlotych galotow"
 632 cc           write(iout,*) "efree=",efree,iii
 633 cc           endif
 634           endif
 635           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 636 #ifdef DEBUG
 637           do k=1,2*nres
 638             do l=1,3
 639               c(l,k)=csingle(l,k)
 640             enddo
 641           enddo
 642           call int_from_cart1(.false.)
 643           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
 644           write (iout,*) "The Cartesian geometry is:"
 645           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 646           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 647           write (iout,*) "The internal geometry is:"
 648           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 649           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 650           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 651           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 652           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 653           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 654           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 655           write (iout,'(f10.5,i5)') rmsdev,iscor
 656 #endif
 657         enddo ! i
 658         write (iout,*) iii," conformations (from",indstart(j)," to",
 659      &   indend(j),") read from ",
 660      &   bxname(:ilen(bxname))
 661         close (ientin,status="delete")
 662       enddo ! j
 663       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
 664 #if (defined(AIX) && !defined(JUBL))
 665       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
 666 #else
 667       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
 668 #endif
 669 #endif
 670       return
 671   101 write (iout,*) "Error in scratchfile."
 672       call flush(iout)
 673       return1
 674       end
 675 c-------------------------------------------------------------------------------
 676       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 677       implicit none
 678       include "DIMENSIONS"
 679       include "DIMENSIONS.ZSCOPT"
 680       include "DIMENSIONS.FREE"
 681       include "DIMENSIONS.COMPAR"
 682 #ifdef MPI
 683       include "mpif.h"
 684       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 685       include "COMMON.MPI"
 686 #endif
 687       include "COMMON.CONTROL"
 688       include "COMMON.CHAIN"
 689       include "COMMON.IOUNITS"
 690       include "COMMON.PROTFILES"
 691       include "COMMON.NAMES"
 692       include "COMMON.VAR"
 693       include "COMMON.SBRIDGE"
 694       include "COMMON.GEO"
 695       include "COMMON.FFIELD"
 696       include "COMMON.ENEPS"
 697       include "COMMON.LOCAL"
 698       include "COMMON.WEIGHTS"
 699       include "COMMON.INTERACT"
 700       include "COMMON.FREE"
 701       include "COMMON.ENERGIES"
 702       include "COMMON.COMPAR"
 703       include "COMMON.PROT"
 704       integer i,j,itmp,iscor,iret,ixdrf
 705       double precision rmsdev,efree,eini
 706       real*4 csingle(3,maxres2),xoord(3,maxres2+2)
 707       real*4 prec
 708
 709 c      write (iout,*) "cxwrite"
 710 c      call flush(iout)
 711       prec=10000.0
 712       do i=1,nres
 713        do j=1,3
 714         xoord(j,i)=csingle(j,i)
 715        enddo
 716       enddo
 717       do i=nnt,nct
 718        do j=1,3
 719         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
 720        enddo
 721       enddo
 722
 723       itmp=nres+nct-nnt+1
 724
 725 c      write (iout,*) "itmp",itmp
 726 c      call flush(iout)
 727 c       write (iout,*) "CNZ",eini,dyn_ss
 728 #if (defined(AIX) && !defined(JUBL))
 729       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 730
 731 c      write (iout,*) "xdrf3dfcoord"
 732 c      call flush(iout)
 733       call xdrfint_(ixdrf, nss, iret)
 734       do j=1,nss
 735 cc       if (dyn_ss) then
 736 cc       call xdrfint_(ixdrf, idssb(j)+nres, iret)
 737 cc        call xdrfint_(ixdrf, jdssb(j)+nres, iret)
 738 cc       else
 739         call xdrfint_(ixdrf, ihpb(j), iret)
 740         call xdrfint_(ixdrf, jhpb(j), iret)
 741 cc       endif
 742       enddo
 743       call xdrffloat_(ixdrf,real(eini),iret)
 744       call xdrffloat_(ixdrf,real(efree),iret)
 745       call xdrffloat_(ixdrf,real(rmsdev),iret)
 746       call xdrfint_(ixdrf,iscor,iret)
 747 #else
 748       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 749
 750       call xdrfint(ixdrf, nss, iret)
 751       do j=1,nss
 752 cc       if (dyn_ss) then
 753 cc        call xdrfint(ixdrf, idssb(j), iret)
 754 cc        call xdrfint(ixdrf, jdssb(j), iret)
 755 cc        idssb(j)=idssb(j)-nres
 756 cc       jdssb(j)=jdssb(j)-nres
 757 cc       else
 758         call xdrfint(ixdrf, ihpb(j), iret)
 759         call xdrfint(ixdrf, jhpb(j), iret)
 760 cc       endif
 761       enddo
 762       call xdrffloat(ixdrf,real(eini),iret)
 763       call xdrffloat(ixdrf,real(efree),iret)
 764       call xdrffloat(ixdrf,real(rmsdev),iret)
 765       call xdrfint(ixdrf,iscor,iret)
 766 #endif
 767
 768       return
 769       end
 770 c------------------------------------------------------------------------------
 771       logical function conf_check(ii,iprint)
 772       implicit none
 773       include "DIMENSIONS"
 774       include "DIMENSIONS.ZSCOPT"
 775       include "DIMENSIONS.FREE"
 776 #ifdef MPI
 777       include "mpif.h"
 778       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 779       include "COMMON.MPI"
 780 #endif
 781       include "COMMON.CHAIN"
 782       include "COMMON.IOUNITS"
 783       include "COMMON.PROTFILES"
 784       include "COMMON.NAMES"
 785       include "COMMON.VAR"
 786       include "COMMON.SBRIDGE"
 787       include "COMMON.GEO"
 788       include "COMMON.FFIELD"
 789       include "COMMON.ENEPS"
 790       include "COMMON.LOCAL"
 791       include "COMMON.WEIGHTS"
 792       include "COMMON.INTERACT"
 793       include "COMMON.FREE"
 794       include "COMMON.ENERGIES"
 795       include "COMMON.CONTROL"
 796       include "COMMON.TORCNSTR"
 797       integer j,k,l,ii,itj,iprint
 798       if (.not.check_conf) then
 799         conf_check=.true.
 800         return
 801       endif
 802       call int_from_cart1(.false.)
 803       do j=nnt+1,nct
 804         if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
 805           if (iprint.gt.0)
 806      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
 807      &      " for conformation",ii
 808           if (iprint.gt.1) then
 809             write (iout,*) "The Cartesian geometry is:"
 810             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 811             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 812             write (iout,*) "The internal geometry is:"
 813             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 814             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 815             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 816             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 817             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 818             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 819           endif
 820           if (iprint.gt.0) write (iout,*)
 821      &      "This conformation WILL NOT be added to the database."
 822           conf_check=.false.
 823           return
 824         endif
 825       enddo
 826       do j=nnt,nct
 827         itj=itype(j)
 828         if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
 829           if (iprint.gt.0)
 830      &    write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),
 831      &     " for conformation",ii
 832           if (iprint.gt.1) then
 833             write (iout,*) "The Cartesian geometry is:"
 834             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 835             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 836             write (iout,*) "The internal geometry is:"
 837             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 838             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 839             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 840             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 841             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 842             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 843           endif
 844           if (iprint.gt.0) write (iout,*)
 845      &      "This conformation WILL NOT be added to the database."
 846           conf_check=.false.
 847           return
 848         endif
 849       enddo
 850       do j=3,nres
 851         if (theta(j).le.0.0d0) then
 852           if (iprint.gt.0)
 853      &    write (iout,*) "Zero theta angle(s) in conformation",ii
 854           if (iprint.gt.1) then
 855             write (iout,*) "The Cartesian geometry is:"
 856             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 857             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 858             write (iout,*) "The internal geometry is:"
 859             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 860             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 861             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 862             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 863             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 864             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 865           endif
 866           if (iprint.gt.0) write (iout,*)
 867      &      "This conformation WILL NOT be added to the database."
 868           conf_check=.false.
 869           return
 870         endif
 871         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 872       enddo
 873       conf_check=.true.
 874 c      write (iout,*) "conf_check passed",ii
 875       return
 876       end