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