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