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