source/wham/src/enecalc1.F

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