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