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