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