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