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