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