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