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