source/unres/src-5hdiag-tmp/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.SPLITELE'
  12       logical file_exist
  13 C Read job setup parameters
  14       call read_control
  15 C Read force-field parameters except weights
  16       call parmread
  17 C Read control parameters for energy minimzation if required
  18       if (minim) call read_minim
  19 C Read MCM control parameters if required
  20       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  21 C Read MD control parameters if reqjuired
  22       if (modecalc.eq.12) call read_MDpar
  23 C Read MREMD control parameters if required
  24       if (modecalc.eq.14) then
  25          call read_MDpar
  26          call read_REMDpar
  27       endif
  28 C Read MUCA control parameters if required
  29       if (lmuca) call read_muca
  30 C Read CSA control parameters if required (from fort.40 if exists
  31 C otherwise from general input file)
  32       if (modecalc.eq.8) then
  33        inquire (file="fort.40",exist=file_exist)
  34        if (.not.file_exist) call csaread
  35       endif
  36 cfmc      if (modecalc.eq.10) call mcmfread
  37 c Read energy-term weights and disulfide parameters
  38       call weightread
  39 C Read molecule information, molecule geometry, energy-term weights, and
  40 C restraints if requested
  41       call molread
  42 C Print restraint information
  43 #ifdef MPI
  44       if (.not. out1file .or. me.eq.king) then
  45 #endif
  46       if (nhpb.gt.nss)  then
  47         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  48      &  "The following",nhpb-nss,
  49      &  " distance restraints have been imposed:",
  50      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  51      &  "  score"," type"
  52         do i=nss+1,nhpb
  53           write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
  54      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
  55      &  irestr_type(i)
  56         enddo
  57       endif
  58       if (npeak.gt.0) then
  59         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  60      &  "The following",npeak,
  61      &  " NMR peak restraints have been imposed:",
  62      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  63      &  "  score"," type"," ipeak"
  64         do i=1,npeak
  65           do j=ipeak(1,i),ipeak(2,i)
  66             write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
  67      &       jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
  68      &       forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
  69           enddo
  70         enddo
  71         write (iout,*) "The ipeak array"
  72         do i=1,npeak
  73           write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
  74         enddo
  75       endif
  76 #ifdef MPI
  77       endif
  78 #endif
  79 c      print *,"Processor",myrank," leaves READRTNS"
  80       return
  81       end
  82 C-------------------------------------------------------------------------------
  83       subroutine read_control
  84 C
  85 C Read contorl data
  86 C
  87       implicit real*8 (a-h,o-z)
  88       include 'DIMENSIONS'
  89 #ifdef MP
  90       include 'mpif.h'
  91       logical OKRandom, prng_restart
  92       real*8  r1
  93 #endif
  94       include 'COMMON.IOUNITS'
  95       include 'COMMON.TIME1'
  96       include 'COMMON.THREAD'
  97       include 'COMMON.SBRIDGE'
  98       include 'COMMON.CONTROL'
  99       include 'COMMON.MCM'
 100       include 'COMMON.MAP'
 101       include 'COMMON.HEADER'
 102       include 'COMMON.CSA'
 103       include 'COMMON.CHAIN'
 104       include 'COMMON.MUCA'
 105       include 'COMMON.MD'
 106       include 'COMMON.FFIELD'
 107       include 'COMMON.INTERACT'
 108       include 'COMMON.SETUP'
 109       include 'COMMON.SPLITELE'
 110       include 'COMMON.SHIELD'
 111       include 'COMMON.GEO'
 112       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
 113       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
 114       character*80 ucase
 115       character*320 controlcard
 116
 117       nglob_csa=0
 118       eglob_csa=1d99
 119       nmin_csa=0
 120       read (INP,'(a)') titel
 121       call card_concat(controlcard)
 122 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
 123 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 124       call reada(controlcard,'SEED',seed,0.0D0)
 125       call random_init(seed)
 126 C Set up the time limit (caution! The time must be input in minutes!)
 127       read_cart=index(controlcard,'READ_CART').gt.0
 128       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 129 C this variable with_theta_constr is the variable which allow to read and execute the
 130 C constrains on theta angles WITH_THETA_CONSTR is the keyword
 131       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 132       write (iout,*) "with_theta_constr ",with_theta_constr
 133       call readi(controlcard,'NSAXS',nsaxs,0)
 134       call readi(controlcard,'SAXS_MODE',saxs_mode,0)
 135       call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
 136       call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
 137       write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
 138      &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
 139       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
 140       read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
 141       out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
 142       out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
 143       call readi(controlcard,'SYM',symetr,1)
 144       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
 145       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 146       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 147       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 148       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 149       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 150       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 151       call reada(controlcard,'DRMS',drms,0.1D0)
 152       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 153        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 154        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 155        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 156        write (iout,'(a,f10.1)')'DRMS    = ',drms
 157        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 158        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 159       endif
 160       call readi(controlcard,'NZ_START',nz_start,0)
 161       call readi(controlcard,'NZ_END',nz_end,0)
 162 c      call readi(controlcard,'IZ_SC',iz_sc,0)
 163       timlim=60.0D0*timlim
 164       safety = 60.0d0*safety
 165       timem=timlim
 166       modecalc=0
 167       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 168       minim=(index(controlcard,'MINIMIZE').gt.0)
 169       dccart=(index(controlcard,'CART').gt.0)
 170       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 171       overlapsc=.not.overlapsc
 172       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 173       searchsc=.not.searchsc
 174       sideadd=(index(controlcard,'SIDEADD').gt.0)
 175       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 176       mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
 177       outpdb=(index(controlcard,'PDBOUT').gt.0)
 178       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 179       pdbref=(index(controlcard,'PDBREF').gt.0)
 180       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 181       indpdb=index(controlcard,'PDBSTART')
 182       extconf=(index(controlcard,'EXTCONF').gt.0)
 183       AFMlog=(index(controlcard,'AFM'))
 184       selfguide=(index(controlcard,'SELFGUIDE'))
 185 c      print *,'AFMlog',AFMlog,selfguide,"KUPA"
 186       call readi(controlcard,'TUBEMOD',tubelog,0)
 187 c      write (iout,*) TUBElog,"TUBEMODE"
 188       call readi(controlcard,'IPRINT',iprint,0)
 189 C SHIELD keyword sets if the shielding effect of side-chains is used
 190 C 0 denots no shielding is used all peptide are equally despite the
 191 C solvent accesible area
 192 C 1 the newly introduced function
 193 C 2 reseved for further possible developement
 194       call readi(controlcard,'SHIELD',shield_mode,0)
 195 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 196         write(iout,*) "shield_mode",shield_mode
 197 C      endif
 198       call readi(controlcard,'TORMODE',tor_mode,0)
 199 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 200         write(iout,*) "torsional and valence angle mode",tor_mode
 201       call readi(controlcard,'MAXGEN',maxgen,10000)
 202       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 203       call readi(controlcard,"KDIAG",kdiag,0)
 204       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 205       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 206      & write (iout,*) "RESCALE_MODE",rescale_mode
 207       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 208       if (index(controlcard,'REGULAR').gt.0.0D0) then
 209         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 210         modecalc=1
 211         refstr=.true.
 212       endif
 213       if (refstr) then
 214       if (index(controlcard,"CASC").gt.0) then
 215         iz_sc=1
 216 c      else if (index(controlcard,"CAONLY").gt.0) then
 217 c        iz_sc=0
 218       else if (index(controlcard,"SCONLY").gt.0) then
 219         iz_sc=2
 220       else
 221         iz_sc=0
 222 c        write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
 223 c     &    " specifying REFSTR, PDBREF or REGULAR."
 224 c        stop
 225       endif
 226       endif
 227       if (index(controlcard,'CHECKGRAD').gt.0) then
 228         modecalc=5
 229         if (index(controlcard,' CART').gt.0) then
 230           icheckgrad=1
 231         elseif (index(controlcard,'CARINT').gt.0) then
 232           icheckgrad=2
 233         else
 234           icheckgrad=3
 235         endif
 236         call reada(controlcard,'DELTA',aincr,1.0d-7)
 237 c        write (iout,*) "icheckgrad",icheckgrad
 238       elseif (index(controlcard,'THREAD').gt.0) then
 239         modecalc=2
 240         call readi(controlcard,'THREAD',nthread,0)
 241         if (nthread.gt.0) then
 242           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 243         else
 244           if (fg_rank.eq.0)
 245      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 246           stop 'Error termination in Read_Control.'
 247         endif
 248       else if (index(controlcard,'MCMA').gt.0) then
 249         modecalc=3
 250       else if (index(controlcard,'MCEE').gt.0) then
 251         modecalc=6
 252       else if (index(controlcard,'MULTCONF').gt.0) then
 253         modecalc=4
 254       else if (index(controlcard,'MAP').gt.0) then
 255         modecalc=7
 256         call readi(controlcard,'MAP',nmap,0)
 257       else if (index(controlcard,'CSA').gt.0) then
 258         modecalc=8
 259 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 260 crc
 261 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 262 crc
 263 crc        modecalc=9
 264 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 265 cfmc        modecalc=10
 266       else if (index(controlcard,'SOFTREG').gt.0) then
 267         modecalc=11
 268       else if (index(controlcard,'CHECK_BOND').gt.0) then
 269         modecalc=-1
 270       else if (index(controlcard,'TEST').gt.0) then
 271         modecalc=-2
 272       else if (index(controlcard,'MD').gt.0) then
 273         modecalc=12
 274       else if (index(controlcard,'RE ').gt.0) then
 275         modecalc=14
 276       endif
 277
 278       lmuca=index(controlcard,'MUCA').gt.0
 279       call readi(controlcard,'MUCADYN',mucadyn,0)
 280       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 281       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 282      & then
 283        write (iout,*) 'MUCADYN=',mucadyn
 284        write (iout,*) 'MUCASMOOTH=',muca_smooth
 285       endif
 286
 287       iscode=index(controlcard,'ONE_LETTER')
 288       indphi=index(controlcard,'PHI')
 289       indback=index(controlcard,'BACK')
 290       iranconf=index(controlcard,'RAND_CONF')
 291       i2ndstr=index(controlcard,'USE_SEC_PRED')
 292       gradout=index(controlcard,'GRADOUT').gt.0
 293       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 294 C  DISTCHAINMAX become obsolete for periodic boundry condition
 295       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 296 C Reading the dimensions of box in x,y,z coordinates
 297       call reada(controlcard,'BOXX',boxxsize,100.0d0)
 298       call reada(controlcard,'BOXY',boxysize,100.0d0)
 299       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 300 c Cutoff range for interactions
 301       call reada(controlcard,"R_CUT",r_cut,15.0d0)
 302       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 303       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 304       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 305       if (lipthick.gt.0.0d0) then
 306        bordliptop=(boxzsize+lipthick)/2.0
 307        bordlipbot=bordliptop-lipthick
 308 C      endif
 309       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 310      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 311       buflipbot=bordlipbot+lipbufthick
 312       bufliptop=bordliptop-lipbufthick
 313       if ((lipbufthick*2.0d0).gt.lipthick)
 314      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 315       endif
 316       write(iout,*) "bordliptop=",bordliptop
 317       write(iout,*) "bordlipbot=",bordlipbot
 318       write(iout,*) "bufliptop=",bufliptop
 319       write(iout,*) "buflipbot=",buflipbot
 320       write (iout,*) "SHIELD MODE",shield_mode
 321       if (TUBElog.gt.0) then
 322        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
 323        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
 324        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
 325        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
 326        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
 327        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
 328        buftubebot=bordtubebot+tubebufthick
 329        buftubetop=bordtubetop-tubebufthick
 330       endif
 331 c      if (shield_mode.gt.0) then
 332 c      pi=3.141592d0
 333 C VSolvSphere the volume of solving sphere
 334 C      print *,pi,"pi"
 335 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
 336 C there will be no distinction between proline peptide group and normal peptide
 337 C group in case of shielding parameters
 338 c      write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
 339 c      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
 340 c      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
 341 c      write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
 342 c     &  VSolvSphere_div
 343 C long axis of side chain
 344 c      do i=1,ntyp
 345 c      long_r_sidechain(i)=vbldsc0(1,i)
 346 c      short_r_sidechain(i)=sigma0(i)
 347 c      enddo
 348 c      buff_shield=1.0d0
 349 c      endif
 350       if (me.eq.king .or. .not.out1file )
 351      &  write (iout,*) "DISTCHAINMAX",distchainmax
 352
 353       if(me.eq.king.or..not.out1file)
 354      & write (iout,'(2a)') diagmeth(kdiag),
 355      &  ' routine used to diagonalize matrices.'
 356       return
 357       end
 358 c--------------------------------------------------------------------------
 359       subroutine read_REMDpar
 360 C
 361 C Read REMD settings
 362 C
 363       implicit real*8 (a-h,o-z)
 364       include 'DIMENSIONS'
 365       include 'COMMON.IOUNITS'
 366       include 'COMMON.TIME1'
 367       include 'COMMON.MD'
 368 #ifndef LANG0
 369       include 'COMMON.LANGEVIN'
 370 #else
 371       include 'COMMON.LANGEVIN.lang0'
 372 #endif
 373       include 'COMMON.INTERACT'
 374       include 'COMMON.NAMES'
 375       include 'COMMON.GEO'
 376       include 'COMMON.REMD'
 377       include 'COMMON.CONTROL'
 378       include 'COMMON.SETUP'
 379       character*80 ucase
 380       character*320 controlcard
 381       character*3200 controlcard1
 382       integer iremd_m_total
 383
 384       if(me.eq.king.or..not.out1file)
 385      & write (iout,*) "REMD setup"
 386
 387       call card_concat(controlcard)
 388       call readi(controlcard,"NREP",nrep,3)
 389       call readi(controlcard,"NSTEX",nstex,1000)
 390       call reada(controlcard,"RETMIN",retmin,10.0d0)
 391       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 392       mremdsync=(index(controlcard,'SYNC').gt.0)
 393       call readi(controlcard,"NSYN",i_sync_step,100)
 394       restart1file=(index(controlcard,'REST1FILE').gt.0)
 395       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 396       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 397       if(max_cache_traj_use.gt.max_cache_traj)
 398      &           max_cache_traj_use=max_cache_traj
 399       if(me.eq.king.or..not.out1file) then
 400 cd       if (traj1file) then
 401 crc caching is in testing - NTWX is not ignored
 402 cd        write (iout,*) "NTWX value is ignored"
 403 cd        write (iout,*) "  trajectory is stored to one file by master"
 404 cd        write (iout,*) "  before exchange at NSTEX intervals"
 405 cd       endif
 406        write (iout,*) "NREP= ",nrep
 407        write (iout,*) "NSTEX= ",nstex
 408        write (iout,*) "SYNC= ",mremdsync
 409        write (iout,*) "NSYN= ",i_sync_step
 410        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 411       endif
 412       remd_tlist=.false.
 413       if (index(controlcard,'TLIST').gt.0) then
 414          remd_tlist=.true.
 415          call card_concat(controlcard1)
 416          read(controlcard1,*) (remd_t(i),i=1,nrep)
 417          if(me.eq.king.or..not.out1file)
 418      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 419       endif
 420       remd_mlist=.false.
 421       if (index(controlcard,'MLIST').gt.0) then
 422          remd_mlist=.true.
 423          call card_concat(controlcard1)
 424          read(controlcard1,*) (remd_m(i),i=1,nrep)
 425          if(me.eq.king.or..not.out1file) then
 426           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 427           iremd_m_total=0
 428           do i=1,nrep
 429            iremd_m_total=iremd_m_total+remd_m(i)
 430           enddo
 431            write (iout,*) 'Total number of replicas ',iremd_m_total
 432           endif
 433          endif
 434         if (adaptive) then
 435           write (iout,*)
 436      &      "Adaptive (PMF-biased) umbrella sampling will be run"
 437           call PMFread
 438         endif
 439       if(me.eq.king.or..not.out1file)
 440      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 441       return
 442       end
 443 c--------------------------------------------------------------------------
 444       subroutine read_MDpar
 445 C
 446 C Read MD settings
 447 C
 448       implicit real*8 (a-h,o-z)
 449       include 'DIMENSIONS'
 450       include 'COMMON.IOUNITS'
 451       include 'COMMON.TIME1'
 452       include 'COMMON.MD'
 453       include 'COMMON.QRESTR'
 454 #ifndef LANG0
 455       include 'COMMON.LANGEVIN'
 456 #else
 457       include 'COMMON.LANGEVIN.lang0'
 458 #endif
 459       include 'COMMON.INTERACT'
 460       include 'COMMON.NAMES'
 461       include 'COMMON.GEO'
 462       include 'COMMON.SETUP'
 463       include 'COMMON.CONTROL'
 464       include 'COMMON.SPLITELE'
 465       include 'COMMON.FFIELD'
 466       character*80 ucase
 467       character*320 controlcard
 468
 469       call card_concat(controlcard)
 470       call readi(controlcard,"NSTEP",n_timestep,1000000)
 471       call readi(controlcard,"NTWE",ntwe,100)
 472       call readi(controlcard,"NTWX",ntwx,1000)
 473       call reada(controlcard,"DT",d_time,1.0d-1)
 474       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 475       call reada(controlcard,"DAMAX",damax,1.0d1)
 476       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 477       call readi(controlcard,"LANG",lang,0)
 478       RESPA = index(controlcard,"RESPA") .gt. 0
 479       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 480       ntime_split0=ntime_split
 481       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 482       ntime_split0=ntime_split
 483 c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
 484 c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 485       rest = index(controlcard,"REST").gt.0
 486       tbf = index(controlcard,"TBF").gt.0
 487       usampl = index(controlcard,"USAMPL").gt.0
 488       scale_umb =  index(controlcard,"SCALE_UMB").gt.0
 489       adaptive = index(controlcard,"ADAPTIVE").gt.0
 490       mdpdb = index(controlcard,"MDPDB").gt.0
 491       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 492       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 493       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 494       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 495       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 496       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 497       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 498       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 499       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 500       large = index(controlcard,"LARGE").gt.0
 501       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 502       rattle = index(controlcard,"RATTLE").gt.0
 503       preminim = index(controlcard,"PREMINIM").gt.0
 504       if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
 505         preminim=.true.
 506         dccart=.true.
 507       endif
 508       write (iout,*) "PREMINIM ",preminim
 509       if (preminim) then
 510         if (index(controlcard,'CART').gt.0) dccart=.true.
 511         write (iout,*) "dccart ",dccart
 512         write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
 513         if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
 514       endif
 515 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 516       nset=0
 517       if(usampl) then
 518         write (iout,*) "Umbrella sampling will be run"
 519         if (scale_umb.and.adaptive) then
 520           write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
 521           write (iout,*) "Select one of those and re-run the job."
 522           stop
 523         endif
 524         if (scale_umb) write (iout,*)
 525      &"Umbrella-restraint force constants will be scaled by temperature"
 526         call read_fragments
 527       endif
 528
 529       if(me.eq.king.or..not.out1file) then
 530        write (iout,*)
 531        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 532        write (iout,*)
 533        write (iout,'(a)') "The units are:"
 534        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 535        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 536      &  " acceleration: angstrom/(48.9 fs)**2"
 537        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 538        write (iout,*)
 539        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 540        write (iout,'(a60,f10.5,a)')
 541      &  "Initial time step of numerical integration:",d_time,
 542      &  " natural units"
 543        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 544        if (RESPA) then
 545         write (iout,'(2a,i4,a)')
 546      &    "A-MTS algorithm used; initial time step for fast-varying",
 547      &    " short-range forces split into",ntime_split," steps."
 548         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 549      &   r_cut," lambda",rlamb
 550        endif
 551        write (iout,'(2a,f10.5)')
 552      &  "Maximum acceleration threshold to reduce the time step",
 553      &  "/increase split number:",damax
 554        write (iout,'(2a,f10.5)')
 555      &  "Maximum predicted energy drift to reduce the timestep",
 556      &  "/increase split number:",edriftmax
 557        write (iout,'(a60,f10.5)')
 558      & "Maximum velocity threshold to reduce velocities:",dvmax
 559        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 560        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 561        if (rattle) write (iout,'(a60)')
 562      &  "Rattle algorithm used to constrain the virtual bonds"
 563       endif
 564       reset_fricmat=1000
 565       if (lang.gt.0) then
 566         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 567         call reada(controlcard,"RWAT",rwat,1.4d0)
 568         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 569         surfarea=index(controlcard,"SURFAREA").gt.0
 570         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 571         if(me.eq.king.or..not.out1file)then
 572          write (iout,'(/a,$)') "Langevin dynamics calculation"
 573          if (lang.eq.1) then
 574           write (iout,'(a/)')
 575      &      " with direct integration of Langevin equations"
 576          else if (lang.eq.2) then
 577           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 578          else if (lang.eq.3) then
 579           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 580          else if (lang.eq.4) then
 581           write (iout,'(a/)') " in overdamped mode"
 582          else
 583           write (iout,'(//a,i5)')
 584      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 585           stop
 586          endif
 587          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 588          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 589          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 590          write (iout,'(a60,f10.5)')
 591      &   "Scaling factor of the friction forces:",scal_fric
 592          if (surfarea) write (iout,'(2a,i10,a)')
 593      &     "Friction coefficients will be scaled by solvent-accessible",
 594      &     " surface area every",reset_fricmat," steps."
 595         endif
 596 c Calculate friction coefficients and bounds of stochastic forces
 597         eta=6*pi*cPoise*etawat
 598         if(me.eq.king.or..not.out1file)
 599      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 600      &   ,eta
 601         gamp=scal_fric*(pstok+rwat)*eta
 602         stdfp=dsqrt(2*Rb*t_bath/d_time)
 603         do i=1,ntyp
 604           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 605           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 606         enddo
 607         if(me.eq.king.or..not.out1file)then
 608          write (iout,'(/2a/)')
 609      &   "Radii of site types and friction coefficients and std's of",
 610      &   " stochastic forces of fully exposed sites"
 611          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 612          do i=1,ntyp
 613           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 614      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 615          enddo
 616         endif
 617       else if (tbf) then
 618         if(me.eq.king.or..not.out1file)then
 619          write (iout,'(a)') "Berendsen bath calculation"
 620          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 621          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 622          if (reset_moment)
 623      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 624      &   count_reset_moment," steps"
 625          if (reset_vel)
 626      &    write (iout,'(a,i10,a)')
 627      &    "Velocities will be reset at random every",count_reset_vel,
 628      &   " steps"
 629         endif
 630       else
 631         if(me.eq.king.or..not.out1file)
 632      &   write (iout,'(a31)') "Microcanonical mode calculation"
 633       endif
 634       if(me.eq.king.or..not.out1file)then
 635        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 636        if (usampl) then
 637           write(iout,*) "MD running with constraints."
 638           write(iout,*) "Equilibration time ", eq_time, " mtus."
 639           write(iout,*) "Constraining ", nfrag," fragments."
 640           write(iout,*) "Length of each fragment, weight and q0:"
 641           do iset=1,nset
 642            write (iout,*) "Set of restraints #",iset
 643            do i=1,nfrag
 644               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 645      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 646            enddo
 647            write(iout,*) "constraints between ", npair, "fragments."
 648            write(iout,*) "constraint pairs, weights and q0:"
 649            do i=1,npair
 650             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 651      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 652            enddo
 653            write(iout,*) "angle constraints within ", nfrag_back,
 654      &      "backbone fragments."
 655            if (loc_qlike) then
 656            write(iout,*) "fragment, weights, q0:"
 657            do i=1,nfrag_back
 658             write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
 659      &         ifrag_back(2,i,iset),
 660      &         wfrag_back(1,i,iset),qin_back(1,i,iset),
 661      &         wfrag_back(2,i,iset),qin_back(2,i,iset),
 662      &         wfrag_back(3,i,iset),qin_back(3,i,iset)
 663            enddo
 664            else
 665            write(iout,*) "fragment, weights:"
 666            do i=1,nfrag_back
 667             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 668      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 669      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 670            enddo
 671            endif
 672           enddo
 673         iset=mod(kolor,nset)+1
 674        endif
 675       endif
 676       if(me.eq.king.or..not.out1file)
 677      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 678       return
 679       end
 680 c------------------------------------------------------------------------------
 681       subroutine molread
 682 C
 683 C Read molecular data.
 684 C
 685       implicit real*8 (a-h,o-z)
 686       include 'DIMENSIONS'
 687 #ifdef MPI
 688       include 'mpif.h'
 689       integer error_msg
 690 #endif
 691       include 'COMMON.IOUNITS'
 692       include 'COMMON.GEO'
 693       include 'COMMON.VAR'
 694       include 'COMMON.INTERACT'
 695       include 'COMMON.LOCAL'
 696       include 'COMMON.NAMES'
 697       include 'COMMON.CHAIN'
 698       include 'COMMON.FFIELD'
 699       include 'COMMON.SBRIDGE'
 700       include 'COMMON.HEADER'
 701       include 'COMMON.CONTROL'
 702       include 'COMMON.DBASE'
 703       include 'COMMON.THREAD'
 704       include 'COMMON.CONTACTS'
 705       include 'COMMON.TORCNSTR'
 706       include 'COMMON.TIME1'
 707       include 'COMMON.BOUNDS'
 708       include 'COMMON.HOMOLOGY'
 709       include 'COMMON.SETUP'
 710       include 'COMMON.SHIELD'
 711       character*4 sequence(maxres)
 712       integer rescode
 713       double precision x(maxvar)
 714       character*256 pdbfile
 715       character*400 weightcard
 716       character*80 weightcard_t,ucase
 717       dimension itype_pdb(maxres)
 718       common /pizda/ itype_pdb
 719       logical seq_comp,fail
 720       double precision energia(0:n_ene)
 721       double precision secprob(3,maxdih_constr)
 722       integer ilen
 723       external ilen
 724 C
 725 C Read PDB structure if applicable
 726 C
 727       if (indpdb.gt.0 .or. pdbref) then
 728         read(inp,'(a)') pdbfile
 729         if(me.eq.king.or..not.out1file)
 730      &   write (iout,'(2a)') 'PDB data will be read from file ',
 731      &   pdbfile(:ilen(pdbfile))
 732         open(ipdbin,file=pdbfile,status='old',err=33)
 733         goto 34
 734   33    write (iout,'(a)') 'Error opening PDB file.'
 735         stop
 736   34    continue
 737         call readpdb
 738         do i=1,2*nres
 739           do j=1,3
 740             crefjlee(j,i)=c(j,i)
 741           enddo
 742         enddo
 743 #ifdef DEBUG
 744         do i=1,nres
 745           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 746      &      (crefjlee(j,i+nres),j=1,3)
 747         enddo
 748 #endif
 749         if(me.eq.king.or..not.out1file)
 750      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 751      &   ' nstart_sup=',nstart_sup
 752         do i=1,nres
 753           itype_pdb(i)=itype(i)
 754         enddo
 755         close (ipdbin)
 756         nnt=nstart_sup
 757         nct=nstart_sup+nsup-1
 758         call contact(.false.,ncont_ref,icont_ref,co)
 759
 760         if (sideadd) then
 761          if(me.eq.king.or..not.out1file)
 762      &    write(iout,*)'Adding sidechains'
 763          maxsi=1000
 764          do i=2,nres-1
 765           iti=itype(i)
 766           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 767             nsi=0
 768             fail=.true.
 769             do while (fail.and.nsi.le.maxsi)
 770               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 771               nsi=nsi+1
 772             enddo
 773             if(fail) write(iout,*)'Adding sidechain failed for res ',
 774      &              i,' after ',nsi,' trials'
 775           endif
 776          enddo
 777         endif
 778       endif
 779       if (indpdb.eq.0) then
 780 C Read sequence if not taken from the pdb file.
 781         read (inp,*) nres
 782 c        print *,'nres=',nres
 783         if (iscode.gt.0) then
 784           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 785         else
 786           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 787         endif
 788 C Convert sequence to numeric code
 789         do i=1,nres
 790           itype(i)=rescode(i,sequence(i),iscode)
 791         enddo
 792 C Assign initial virtual bond lengths
 793         do i=2,nres
 794           vbld(i)=vbl
 795           vbld_inv(i)=vblinv
 796         enddo
 797         do i=2,nres-1
 798           vbld(i+nres)=dsc(iabs(itype(i)))
 799           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 800 c          write (iout,*) "i",i," itype",itype(i),
 801 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 802         enddo
 803       endif
 804 c      print *,nres
 805 c      print '(20i4)',(itype(i),i=1,nres)
 806       do i=1,nres
 807 #ifdef PROCOR
 808         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 809 #else
 810         if (itype(i).eq.ntyp1) then
 811 #endif
 812           itel(i)=0
 813 #ifdef PROCOR
 814         else if (iabs(itype(i+1)).ne.20) then
 815 #else
 816         else if (iabs(itype(i)).ne.20) then
 817 #endif
 818           itel(i)=1
 819         else
 820           itel(i)=2
 821         endif
 822       enddo
 823       if(me.eq.king.or..not.out1file)then
 824        write (iout,*) "ITEL"
 825        do i=1,nres-1
 826          write (iout,*) i,itype(i),itel(i)
 827        enddo
 828        print *,'Call Read_Bridge.'
 829       endif
 830       nnt=1
 831       nct=nres
 832 cd      print *,'NNT=',NNT,' NCT=',NCT
 833       call seq2chains(nres,itype,nchain,chain_length,chain_border,
 834      &  ireschain)
 835       write(iout,*) "nres",nres," nchain",nchain
 836       do i=1,nchain
 837         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
 838      &    chain_border(2,i)
 839       enddo
 840       call chain_symmetry(nchain,nres,itype,chain_border,
 841      &    chain_length,npermchain,tabpermchain)
 842       do i=1,nres
 843         write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
 844      &    ii=1,npermchain)
 845       enddo
 846       write(iout,*) "residue permutations"
 847       do i=1,nres
 848         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
 849       enddo
 850       if (itype(1).eq.ntyp1) nnt=2
 851       if (itype(nres).eq.ntyp1) nct=nct-1
 852 #ifdef DFA
 853       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 854      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 855        call init_dfa_vars
 856        print*, 'init_dfa_vars finished!'
 857        call read_dfa_info
 858        print*, 'read_dfa_info finished!'
 859       endif
 860 #endif
 861       if (pdbref) then
 862         if(me.eq.king.or..not.out1file)
 863      &   write (iout,'(a,i3)') 'nsup=',nsup
 864         nstart_seq=nnt
 865         if (nsup.le.(nct-nnt+1)) then
 866           do i=0,nct-nnt+1-nsup
 867             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 868               nstart_seq=nnt+i
 869               goto 111
 870             endif
 871           enddo
 872           write (iout,'(a)')
 873      &            'Error - sequences to be superposed do not match.'
 874           stop
 875         else
 876           do i=0,nsup-(nct-nnt+1)
 877             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 878      &      then
 879               nstart_sup=nstart_sup+i
 880               nsup=nct-nnt+1
 881               goto 111
 882             endif
 883           enddo
 884           write (iout,'(a)')
 885      &            'Error - sequences to be superposed do not match.'
 886         endif
 887   111   continue
 888         if (nsup.eq.0) nsup=nct-nnt
 889         if (nstart_sup.eq.0) nstart_sup=nnt
 890         if (nstart_seq.eq.0) nstart_seq=nnt
 891         if(me.eq.king.or..not.out1file)
 892      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 893      &                 ' nstart_seq=',nstart_seq
 894       endif
 895       call read_bridge
 896 C 8/13/98 Set limits to generating the dihedral angles
 897       do i=1,nres
 898         phibound(1,i)=-pi
 899         phibound(2,i)=pi
 900       enddo
 901       read (inp,*) ndih_constr
 902       write (iout,*) "ndih_constr",ndih_constr
 903       if (ndih_constr.gt.0) then
 904         raw_psipred=.false.
 905 C        read (inp,*) ftors
 906         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 907      &  i=1,ndih_constr)
 908 #ifdef MPI
 909         if(me.eq.king.or..not.out1file)then
 910 #endif
 911          write (iout,*)
 912      &   'There are',ndih_constr,' restraints on gamma angles.'
 913          do i=1,ndih_constr
 914           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 915      &    ftors(i)
 916          enddo
 917 #ifdef MPI
 918         endif
 919 #endif
 920         do i=1,ndih_constr
 921           phi0(i)=deg2rad*phi0(i)
 922           drange(i)=deg2rad*drange(i)
 923         enddo
 924 C        if(me.eq.king.or..not.out1file)
 925 C     &   write (iout,*) 'FTORS',ftors
 926         do i=1,ndih_constr
 927           ii = idih_constr(i)
 928           phibound(1,ii) = phi0(i)-drange(i)
 929           phibound(2,ii) = phi0(i)+drange(i)
 930         enddo
 931 #ifdef MPI
 932         if (me.eq.king .or. .not.out1file) then
 933 #endif
 934          write (iout,'(a)') 'Boundaries in gamma angle sampling:'
 935          do i=1,nres
 936            write (iout,'(a3,i5,2f10.1)')
 937      &    restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 938          enddo
 939 #ifdef MP
 940         endif
 941 #endif
 942       else if (ndih_constr.lt.0) then
 943         raw_psipred=.true.
 944         call card_concat(weightcard)
 945         call reada(weightcard,"PHIHEL",phihel,50.0D0)
 946         call reada(weightcard,"PHIBET",phibet,180.0D0)
 947         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
 948         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
 949         call reada(weightcard,"WDIHC",wdihc,0.591D0)
 950         write (iout,*) "Weight of dihedral angle restraints",wdihc
 951         read(inp,'(9x,3f7.3)')
 952      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
 953         write (iout,*) "The secprob array"
 954         do i=nnt,nct
 955           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
 956         enddo
 957         ndih_constr=0
 958         do i=nnt+3,nct
 959           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
 960      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
 961             ndih_constr=ndih_constr+1
 962             idih_constr(ndih_constr)=i
 963             sumv=0.0d0
 964             do j=1,3
 965               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
 966               sumv=sumv+vpsipred(j,ndih_constr)
 967             enddo
 968             do j=1,3
 969               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
 970             enddo
 971             phibound(1,ndih_constr)=phihel*deg2rad
 972             phibound(2,ndih_constr)=phibet*deg2rad
 973             sdihed(1,ndih_constr)=sigmahel*deg2rad
 974             sdihed(2,ndih_constr)=sigmabet*deg2rad
 975           endif
 976         enddo
 977 #ifdef MPI
 978         if(me.eq.king.or..not.out1file)then
 979 #endif
 980          write (iout,*)
 981      &   'There are',ndih_constr,
 982      &   ' bimodal restraints on gamma angles.'
 983          do i=1,ndih_constr
 984            write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
 985      &      restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
 986      &      restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
 987      &      phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
 988      &      phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
 989      &      (vpsipred(j,i),j=1,3)
 990          enddo
 991 #ifdef MPI
 992         endif
 993 #endif
 994 c Raw psipred input
 995       endif
 996 C first setting the theta boundaries to 0 to pi
 997 C this mean that there is no energy penalty for any angle occuring this can be applied
 998 C for generate random conformation but is not implemented in this way
 999 C      do i=1,nres
1000 C        thetabound(1,i)=0
1001 C        thetabound(2,i)=pi
1002 C      enddo
1003 C begin reading theta constrains this is quartic constrains allowing to
1004 C have smooth second derivative
1005       if (with_theta_constr) then
1006 C with_theta_constr is keyword allowing for occurance of theta constrains
1007       read (inp,*) ntheta_constr
1008 C ntheta_constr is the number of theta constrains
1009       if (ntheta_constr.gt.0) then
1010 C        read (inp,*) ftors
1011         read (inp,*) (itheta_constr(i),theta_constr0(i),
1012      &  theta_drange(i),for_thet_constr(i),
1013      &  i=1,ntheta_constr)
1014 C the above code reads from 1 to ntheta_constr
1015 C itheta_constr(i) residue i for which is theta_constr
1016 C theta_constr0 the global minimum value
1017 C theta_drange is range for which there is no energy penalty
1018 C for_thet_constr is the force constant for quartic energy penalty
1019 C E=k*x**4
1020         if(me.eq.king.or..not.out1file)then
1021          write (iout,*)
1022      &   'There are',ntheta_constr,' constraints on phi angles.'
1023          do i=1,ntheta_constr
1024           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1025      &    theta_drange(i),
1026      &    for_thet_constr(i)
1027          enddo
1028         endif
1029         do i=1,ntheta_constr
1030           theta_constr0(i)=deg2rad*theta_constr0(i)
1031           theta_drange(i)=deg2rad*theta_drange(i)
1032         enddo
1033 C        if(me.eq.king.or..not.out1file)
1034 C     &   write (iout,*) 'FTORS',ftors
1035 C        do i=1,ntheta_constr
1036 C          ii = itheta_constr(i)
1037 C          thetabound(1,ii) = phi0(i)-drange(i)
1038 C          thetabound(2,ii) = phi0(i)+drange(i)
1039 C        enddo
1040       endif ! ntheta_constr.gt.0
1041       endif! with_theta_constr
1042 C
1043 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1044 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
1045 c--- Zscore rms -------
1046       if (nz_start.eq.0) nz_start=nnt
1047       if (nz_end.eq.0 .and. nsup.gt.0) then
1048         nz_end=nnt+nsup-1
1049       else if (nz_end.eq.0) then
1050         nz_end=nct
1051       endif
1052       if(me.eq.king.or..not.out1file)then
1053        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1054        write (iout,*) 'IZ_SC=',iz_sc
1055       endif
1056 c----------------------
1057       call init_int_table
1058       if (refstr) then
1059         if (.not.pdbref) then
1060           call read_angles(inp,*38)
1061           goto 39
1062    38     write (iout,'(a)') 'Error reading reference structure.'
1063 #ifdef MPI
1064           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1065           stop 'Error reading reference structure'
1066 #endif
1067    39     call chainbuild_extconf
1068           call setup_var
1069 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1070           nstart_sup=nnt
1071           nstart_seq=nnt
1072           nsup=nct-nnt+1
1073           do i=1,2*nres
1074             do j=1,3
1075               cref(j,i)=c(j,i)
1076             enddo
1077           enddo
1078           call contact(.true.,ncont_ref,icont_ref,co)
1079         endif
1080       endif
1081 c        print *, "A TU"
1082       write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1083       call flush(iout)
1084       if (constr_dist.gt.0) call read_dist_constr
1085       write (iout,*) "After read_dist_constr nhpb",nhpb
1086       if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1087       call hpb_partition
1088       call NMRpeak_partition
1089       if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1090       if (pdbref) then
1091         if(me.eq.king.or..not.out1file)
1092      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1093         do i=1,ncont_ref
1094           do j=1,2
1095             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1096           enddo
1097           if(me.eq.king.or..not.out1file)
1098      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1099      &     icont_ref(1,i),' ',
1100      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1101         enddo
1102       endif
1103       write (iout,*) "calling read_saxs_consrtr",nsaxs
1104       if (nsaxs.gt.0) call read_saxs_constr
1105
1106       if (constr_homology.gt.0) then
1107         call read_constr_homology
1108         if (indpdb.gt.0 .or. pdbref) then
1109           do i=1,2*nres
1110             do j=1,3
1111               c(j,i)=crefjlee(j,i)
1112               cref(j,i)=crefjlee(j,i)
1113             enddo
1114           enddo
1115         endif
1116 #ifdef DEBUG
1117         write (iout,*) "sc_loc_geom: Array C"
1118         do i=1,nres
1119           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1120      &      (c(j,i+nres),j=1,3)
1121         enddo
1122         write (iout,*) "Array Cref"
1123         do i=1,nres
1124           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1125      &      (cref(j,i+nres),j=1,3)
1126         enddo
1127 #endif
1128        call int_from_cart1(.false.)
1129        call sc_loc_geom(.false.)
1130        do i=1,nres
1131          thetaref(i)=theta(i)
1132          phiref(i)=phi(i)
1133        enddo
1134        do i=1,nres-1
1135          do j=1,3
1136            dc(j,i)=c(j,i+1)-c(j,i)
1137            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1138          enddo
1139        enddo
1140        do i=2,nres-1
1141          do j=1,3
1142            dc(j,i+nres)=c(j,i+nres)-c(j,i)
1143            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1144          enddo
1145        enddo
1146       else
1147         homol_nset=0
1148       endif
1149
1150
1151 C      endif
1152       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1153      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1154      &    modecalc.ne.10) then
1155 C If input structure hasn't been supplied from the PDB file read or generate
1156 C initial geometry.
1157         if (iranconf.eq.0 .and. .not. extconf) then
1158           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1159      &     write (iout,'(a)') 'Initial geometry will be read in.'
1160           if (read_cart) then
1161             read(inp,'(8f10.5)',end=36,err=36)
1162      &       ((c(l,k),l=1,3),k=1,nres),
1163      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1164             write (iout,*) "Exit READ_CART"
1165 c            write (iout,'(8f10.5)')
1166 c     &       ((c(l,k),l=1,3),k=1,nres),
1167 c     &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1168             call cartprint
1169             do i=1,nres-1
1170               do j=1,3
1171                 dc(j,i)=c(j,i+1)-c(j,i)
1172 c                dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1173               enddo
1174             enddo
1175             do i=nnt,nct
1176               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1177                 do j=1,3
1178                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1179 c                  dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1180                 enddo
1181               else
1182                 do j=1,3
1183                   dc(j,i+nres)=0.0d0
1184 c                  dc_norm(j,i+nres)=0.0d0
1185                 enddo
1186               endif
1187             enddo
1188             call int_from_cart1(.true.)
1189             write (iout,*) "Finish INT_TO_CART"
1190 c      write (iout,*) "DC"
1191 c      do i=1,nres
1192 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1193 c     &    (dc(j,i+nres),j=1,3)
1194 c      enddo
1195 c            call cartprint
1196 c            call setup_var
1197 c            return
1198           else
1199             call read_angles(inp,*36)
1200             call chainbuild_extconf
1201           endif
1202           goto 37
1203    36     write (iout,'(a)') 'Error reading angle file.'
1204 #ifdef MPI
1205           call mpi_finalize( MPI_COMM_WORLD,IERR )
1206 #endif
1207           stop 'Error reading angle file.'
1208    37     continue
1209         else if (extconf) then
1210          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1211      &    write (iout,'(a)') 'Extended chain initial geometry.'
1212          do i=3,nres
1213           theta(i)=90d0*deg2rad
1214          enddo
1215          do i=4,nres
1216           phi(i)=180d0*deg2rad
1217          enddo
1218          do i=2,nres-1
1219           alph(i)=110d0*deg2rad
1220          enddo
1221          do i=2,nres-1
1222           omeg(i)=-120d0*deg2rad
1223           if (itype(i).le.0) omeg(i)=-omeg(i)
1224          enddo
1225         call chainbuild_extconf
1226         else
1227           if(me.eq.king.or..not.out1file)
1228      &     write (iout,'(a)') 'Random-generated initial geometry.'
1229
1230
1231 #ifdef MPI
1232           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1233      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1234 #endif
1235             do itrial=1,100
1236               itmp=1
1237               call gen_rand_conf(itmp,*30)
1238               goto 40
1239    30         write (iout,*) 'Failed to generate random conformation',
1240      &          ', itrial=',itrial
1241               write (*,*) 'Processor:',me,
1242      &          ' Failed to generate random conformation',
1243      &          ' itrial=',itrial
1244               call intout
1245
1246 #ifdef AIX
1247               call flush_(iout)
1248 #else
1249               call flush(iout)
1250 #endif
1251             enddo
1252             write (iout,'(a,i3,a)') 'Processor:',me,
1253      &        ' error in generating random conformation.'
1254             write (*,'(a,i3,a)') 'Processor:',me,
1255      &        ' error in generating random conformation.'
1256             call flush(iout)
1257 #ifdef MPI
1258             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1259    40       continue
1260           endif
1261 #else
1262           write (*,'(a)')
1263      &      ' error in generating random conformation.'
1264           stop
1265    40     continue
1266 #endif
1267         endif
1268       elseif (modecalc.eq.4) then
1269         read (inp,'(a)') intinname
1270         open (intin,file=intinname,status='old',err=333)
1271         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1272      &  write (iout,'(a)') 'intinname',intinname
1273         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1274         goto 334
1275   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1276 #ifdef MPI
1277         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1278 #endif
1279         stop 'Error opening angle file.'
1280   334   continue
1281
1282       endif
1283 C Generate distance constraints, if the PDB structure is to be regularized.
1284       if (nthread.gt.0) then
1285         call read_threadbase
1286       endif
1287       call setup_var
1288       if (me.eq.king .or. .not. out1file)
1289      & call intout
1290       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1291         write (iout,'(/a,i3,a)')
1292      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1293         write (iout,'(20i4)') (iss(i),i=1,ns)
1294        if (dyn_ss) then
1295           write(iout,*)"Running with dynamic disulfide-bond formation"
1296        else
1297         write (iout,'(/a/)') 'Pre-formed links are:'
1298         do i=1,nss
1299           i1=ihpb(i)-nres
1300           i2=jhpb(i)-nres
1301           it1=itype(i1)
1302           it2=itype(i2)
1303           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1304      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1305      &    ebr,forcon(i)
1306         enddo
1307         write (iout,'(a)')
1308        endif
1309       endif
1310       if (ns.gt.0.and.dyn_ss) then
1311           do i=nss+1,nhpb
1312             ihpb(i-nss)=ihpb(i)
1313             jhpb(i-nss)=jhpb(i)
1314             forcon(i-nss)=forcon(i)
1315             dhpb(i-nss)=dhpb(i)
1316           enddo
1317           nhpb=nhpb-nss
1318           nss=0
1319           call hpb_partition
1320           do i=1,ns
1321             dyn_ss_mask(iss(i))=.true.
1322           enddo
1323       endif
1324 c --- test
1325 c      call cartprint
1326 c      write (iout,*) "DC"
1327 c      do i=1,nres
1328 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1329 c     &    (dc(j,i+nres),j=1,3)
1330 c      enddo
1331       if (i2ndstr.gt.0) call secstrp2dihc
1332 c      call geom_to_var(nvar,x)
1333 c      call etotal(energia(0))
1334 c      call enerprint(energia(0))
1335 c      call briefout(0,etot)
1336 c      stop
1337 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1338 cd    write (iout,'(a)') 'Variable list:'
1339 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1340 #ifdef MPI
1341       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1342      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1343      &  'Processor',myrank,': end reading molecular data.'
1344 #endif
1345 c      print *,"A TU?"
1346       return
1347       end
1348 c--------------------------------------------------------------------------
1349       logical function seq_comp(itypea,itypeb,length)
1350       implicit none
1351       integer length,itypea(length),itypeb(length)
1352       integer i
1353       do i=1,length
1354         if (itypea(i).ne.itypeb(i)) then
1355           seq_comp=.false.
1356           return
1357         endif
1358       enddo
1359       seq_comp=.true.
1360       return
1361       end
1362 c-----------------------------------------------------------------------------
1363       subroutine read_bridge
1364 C Read information about disulfide bridges.
1365       implicit real*8 (a-h,o-z)
1366       include 'DIMENSIONS'
1367 #ifdef MPI
1368       include 'mpif.h'
1369 #endif
1370       include 'COMMON.IOUNITS'
1371       include 'COMMON.GEO'
1372       include 'COMMON.VAR'
1373       include 'COMMON.INTERACT'
1374       include 'COMMON.LOCAL'
1375       include 'COMMON.NAMES'
1376       include 'COMMON.CHAIN'
1377       include 'COMMON.FFIELD'
1378       include 'COMMON.SBRIDGE'
1379       include 'COMMON.HEADER'
1380       include 'COMMON.CONTROL'
1381       include 'COMMON.DBASE'
1382       include 'COMMON.THREAD'
1383       include 'COMMON.TIME1'
1384       include 'COMMON.SETUP'
1385 C Read bridging residues.
1386       read (inp,*) ns,(iss(i),i=1,ns)
1387       print *,'ns=',ns
1388       if(me.eq.king.or..not.out1file)
1389      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1390 C Check whether the specified bridging residues are cystines.
1391       do i=1,ns
1392         if (itype(iss(i)).ne.1) then
1393           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1394      &   'Do you REALLY think that the residue ',
1395      &    restyp(itype(iss(i))),i,
1396      &   ' can form a disulfide bridge?!!!'
1397           write (*,'(2a,i3,a)')
1398      &   'Do you REALLY think that the residue ',
1399      &    restyp(itype(iss(i))),i,
1400      &   ' can form a disulfide bridge?!!!'
1401 #ifdef MPI
1402          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1403          stop
1404 #endif
1405         endif
1406       enddo
1407 C Read preformed bridges.
1408       if (ns.gt.0) then
1409       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1410       if(fg_rank.eq.0)
1411      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1412       if (nss.gt.0) then
1413         nhpb=nss
1414 C Check if the residues involved in bridges are in the specified list of
1415 C bridging residues.
1416         do i=1,nss
1417           do j=1,i-1
1418             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1419      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1420               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1421      &      ' contains residues present in other pairs.'
1422               write (*,'(a,i3,a)') 'Disulfide pair',i,
1423      &      ' contains residues present in other pairs.'
1424 #ifdef MPI
1425               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1426               stop
1427 #endif
1428             endif
1429           enddo
1430           do j=1,ns
1431             if (ihpb(i).eq.iss(j)) goto 10
1432           enddo
1433           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1434    10     continue
1435           do j=1,ns
1436             if (jhpb(i).eq.iss(j)) goto 20
1437           enddo
1438           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1439    20     continue
1440           dhpb(i)=dbr
1441           forcon(i)=fbr
1442         enddo
1443         do i=1,nss
1444           ihpb(i)=ihpb(i)+nres
1445           jhpb(i)=jhpb(i)+nres
1446         enddo
1447       endif
1448       endif
1449       return
1450       end
1451 c----------------------------------------------------------------------------
1452       subroutine read_x(kanal,*)
1453       implicit real*8 (a-h,o-z)
1454       include 'DIMENSIONS'
1455       include 'COMMON.GEO'
1456       include 'COMMON.VAR'
1457       include 'COMMON.CHAIN'
1458       include 'COMMON.IOUNITS'
1459       include 'COMMON.CONTROL'
1460       include 'COMMON.LOCAL'
1461       include 'COMMON.INTERACT'
1462 c Read coordinates from input
1463 c
1464       read(kanal,'(8f10.5)',end=10,err=10)
1465      &  ((c(l,k),l=1,3),k=1,nres),
1466      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1467       do j=1,3
1468         c(j,nres+1)=c(j,1)
1469         c(j,2*nres)=c(j,nres)
1470       enddo
1471       call int_from_cart1(.false.)
1472       do i=1,nres-1
1473         do j=1,3
1474           dc(j,i)=c(j,i+1)-c(j,i)
1475           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1476         enddo
1477       enddo
1478       do i=nnt,nct
1479         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1480           do j=1,3
1481             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1482             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1483           enddo
1484         endif
1485       enddo
1486
1487       return
1488    10 return1
1489       end
1490 c----------------------------------------------------------------------------
1491       subroutine read_threadbase
1492       implicit real*8 (a-h,o-z)
1493       include 'DIMENSIONS'
1494       include 'COMMON.IOUNITS'
1495       include 'COMMON.GEO'
1496       include 'COMMON.VAR'
1497       include 'COMMON.INTERACT'
1498       include 'COMMON.LOCAL'
1499       include 'COMMON.NAMES'
1500       include 'COMMON.CHAIN'
1501       include 'COMMON.FFIELD'
1502       include 'COMMON.SBRIDGE'
1503       include 'COMMON.HEADER'
1504       include 'COMMON.CONTROL'
1505       include 'COMMON.DBASE'
1506       include 'COMMON.THREAD'
1507       include 'COMMON.TIME1'
1508 C Read pattern database for threading.
1509       read (icbase,*) nseq
1510       do i=1,nseq
1511         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1512      &   nres_base(2,i),nres_base(3,i)
1513         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1514      &   nres_base(1,i))
1515 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1516 c    &   nres_base(2,i),nres_base(3,i)
1517 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1518 c    &   nres_base(1,i))
1519       enddo
1520       close (icbase)
1521       if (weidis.eq.0.0D0) weidis=0.1D0
1522       do i=nnt,nct
1523         do j=i+2,nct
1524           nhpb=nhpb+1
1525           ihpb(nhpb)=i
1526           jhpb(nhpb)=j
1527           forcon(nhpb)=weidis
1528         enddo
1529       enddo
1530       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1531       write (iout,'(a,i5)') 'nexcl: ',nexcl
1532       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1533       return
1534       end
1535 c------------------------------------------------------------------------------
1536       subroutine setup_var
1537       implicit real*8 (a-h,o-z)
1538       include 'DIMENSIONS'
1539       include 'COMMON.IOUNITS'
1540       include 'COMMON.GEO'
1541       include 'COMMON.VAR'
1542       include 'COMMON.INTERACT'
1543       include 'COMMON.LOCAL'
1544       include 'COMMON.NAMES'
1545       include 'COMMON.CHAIN'
1546       include 'COMMON.FFIELD'
1547       include 'COMMON.SBRIDGE'
1548       include 'COMMON.HEADER'
1549       include 'COMMON.CONTROL'
1550       include 'COMMON.DBASE'
1551       include 'COMMON.THREAD'
1552       include 'COMMON.TIME1'
1553 C Set up variable list.
1554       ntheta=nres-2
1555       nphi=nres-3
1556       nvar=ntheta+nphi
1557       nside=0
1558       write (iout,*) "SETUP_VAR ialph"
1559       do i=2,nres-1
1560         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1561           nside=nside+1
1562           ialph(i,1)=nvar+nside
1563           ialph(nside,2)=i
1564         endif
1565       enddo
1566       if (indphi.gt.0) then
1567         nvar=nphi
1568       else if (indback.gt.0) then
1569         nvar=nphi+ntheta
1570       else
1571         nvar=nvar+2*nside
1572       endif
1573       write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1574       return
1575       end
1576 c----------------------------------------------------------------------------
1577       subroutine gen_dist_constr
1578 C Generate CA distance constraints.
1579       implicit real*8 (a-h,o-z)
1580       include 'DIMENSIONS'
1581       include 'COMMON.IOUNITS'
1582       include 'COMMON.GEO'
1583       include 'COMMON.VAR'
1584       include 'COMMON.INTERACT'
1585       include 'COMMON.LOCAL'
1586       include 'COMMON.NAMES'
1587       include 'COMMON.CHAIN'
1588       include 'COMMON.FFIELD'
1589       include 'COMMON.SBRIDGE'
1590       include 'COMMON.HEADER'
1591       include 'COMMON.CONTROL'
1592       include 'COMMON.DBASE'
1593       include 'COMMON.THREAD'
1594       include 'COMMON.TIME1'
1595       dimension itype_pdb(maxres)
1596       common /pizda/ itype_pdb
1597       character*2 iden
1598 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1599 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1600 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1601 cd     & ' nsup',nsup
1602       do i=nstart_sup,nstart_sup+nsup-1
1603 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1604 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1605         do j=i+2,nstart_sup+nsup-1
1606           nhpb=nhpb+1
1607           ihpb(nhpb)=i+nstart_seq-nstart_sup
1608           jhpb(nhpb)=j+nstart_seq-nstart_sup
1609           forcon(nhpb)=weidis
1610           dhpb(nhpb)=dist(i,j)
1611         enddo
1612       enddo
1613 cd      write (iout,'(a)') 'Distance constraints:'
1614 cd      do i=nss+1,nhpb
1615 cd        ii=ihpb(i)
1616 cd        jj=jhpb(i)
1617 cd        iden='CA'
1618 cd        if (ii.gt.nres) then
1619 cd          iden='SC'
1620 cd          ii=ii-nres
1621 cd          jj=jj-nres
1622 cd        endif
1623 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1624 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1625 cd     &  dhpb(i),forcon(i)
1626 cd      enddo
1627       return
1628       end
1629 c----------------------------------------------------------------------------
1630       subroutine map_read
1631       implicit real*8 (a-h,o-z)
1632       include 'DIMENSIONS'
1633       include 'COMMON.MAP'
1634       include 'COMMON.IOUNITS'
1635       character*3 angid(4) /'THE','PHI','ALP','OME'/
1636       character*80 mapcard,ucase
1637       do imap=1,nmap
1638         read (inp,'(a)') mapcard
1639         mapcard=ucase(mapcard)
1640         if (index(mapcard,'PHI').gt.0) then
1641           kang(imap)=1
1642         else if (index(mapcard,'THE').gt.0) then
1643           kang(imap)=2
1644         else if (index(mapcard,'ALP').gt.0) then
1645           kang(imap)=3
1646         else if (index(mapcard,'OME').gt.0) then
1647           kang(imap)=4
1648         else
1649           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1650           stop 'Error - illegal variable spec in MAP card.'
1651         endif
1652         call readi (mapcard,'RES1',res1(imap),0)
1653         call readi (mapcard,'RES2',res2(imap),0)
1654         if (res1(imap).eq.0) then
1655           res1(imap)=res2(imap)
1656         else if (res2(imap).eq.0) then
1657           res2(imap)=res1(imap)
1658         endif
1659         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1660           write (iout,'(a)')
1661      &    'Error - illegal definition of variable group in MAP.'
1662           stop 'Error - illegal definition of variable group in MAP.'
1663         endif
1664         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1665         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1666         call readi(mapcard,'NSTEP',nstep(imap),0)
1667         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1668           write (iout,'(a)')
1669      &     'Illegal boundary and/or step size specification in MAP.'
1670           stop 'Illegal boundary and/or step size specification in MAP.'
1671         endif
1672       enddo ! imap
1673       return
1674       end
1675 c----------------------------------------------------------------------------
1676       subroutine csaread
1677       implicit real*8 (a-h,o-z)
1678       include 'DIMENSIONS'
1679       include 'COMMON.IOUNITS'
1680       include 'COMMON.GEO'
1681       include 'COMMON.CSA'
1682       include 'COMMON.BANK'
1683       include 'COMMON.CONTROL'
1684       character*80 ucase
1685       character*620 mcmcard
1686       call card_concat(mcmcard)
1687
1688       call readi(mcmcard,'NCONF',nconf,50)
1689       call readi(mcmcard,'NADD',nadd,0)
1690       call readi(mcmcard,'JSTART',jstart,1)
1691       call readi(mcmcard,'JEND',jend,1)
1692       call readi(mcmcard,'NSTMAX',nstmax,500000)
1693       call readi(mcmcard,'N0',n0,1)
1694       call readi(mcmcard,'N1',n1,6)
1695       call readi(mcmcard,'N2',n2,4)
1696       call readi(mcmcard,'N3',n3,0)
1697       call readi(mcmcard,'N4',n4,0)
1698       call readi(mcmcard,'N5',n5,0)
1699       call readi(mcmcard,'N6',n6,10)
1700       call readi(mcmcard,'N7',n7,0)
1701       call readi(mcmcard,'N8',n8,0)
1702       call readi(mcmcard,'N9',n9,0)
1703       call readi(mcmcard,'N14',n14,0)
1704       call readi(mcmcard,'N15',n15,0)
1705       call readi(mcmcard,'N16',n16,0)
1706       call readi(mcmcard,'N17',n17,0)
1707       call readi(mcmcard,'N18',n18,0)
1708
1709       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1710
1711       call readi(mcmcard,'NDIFF',ndiff,2)
1712       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1713       call readi(mcmcard,'IS1',is1,1)
1714       call readi(mcmcard,'IS2',is2,8)
1715       call readi(mcmcard,'NRAN0',nran0,4)
1716       call readi(mcmcard,'NRAN1',nran1,2)
1717       call readi(mcmcard,'IRR',irr,1)
1718       call readi(mcmcard,'NSEED',nseed,20)
1719       call readi(mcmcard,'NTOTAL',ntotal,10000)
1720       call reada(mcmcard,'CUT1',cut1,2.0d0)
1721       call reada(mcmcard,'CUT2',cut2,5.0d0)
1722       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1723       call readi(mcmcard,'ICMAX',icmax,3)
1724       call readi(mcmcard,'IRESTART',irestart,0)
1725 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1726       ntbankm=0
1727 c!bankt
1728       call reada(mcmcard,'DELE',dele,20.0d0)
1729       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1730       call readi(mcmcard,'IREF',iref,0)
1731       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1732       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1733       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1734       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1735       write (iout,*) "NCONF_IN",nconf_in
1736       return
1737       end
1738 c----------------------------------------------------------------------------
1739 cfmc      subroutine mcmfread
1740 cfmc      implicit real*8 (a-h,o-z)
1741 cfmc      include 'DIMENSIONS'
1742 cfmc      include 'COMMON.MCMF'
1743 cfmc      include 'COMMON.IOUNITS'
1744 cfmc      include 'COMMON.GEO'
1745 cfmc      character*80 ucase
1746 cfmc      character*620 mcmcard
1747 cfmc      call card_concat(mcmcard)
1748 cfmc
1749 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1750 cfmc      write(iout,*)'MAXRANT=',maxrant
1751 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1752 cfmc      write(iout,*)'MAXFAM=',maxfam
1753 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1754 cfmc      write(iout,*)'NNET1=',nnet1
1755 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1756 cfmc      write(iout,*)'NNET2=',nnet2
1757 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1758 cfmc      write(iout,*)'NNET3=',nnet3
1759 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1760 cfmc      write(iout,*)'ILASTT=',ilastt
1761 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1762 cfmc      write(iout,*)'MAXSTR=',maxstr
1763 cfmc      maxstr_f=maxstr/maxfam
1764 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1765 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1766 cfmc      write(iout,*)'NMCMF=',nmcmf
1767 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1768 cfmc      write(iout,*)'IFOCUS=',ifocus
1769 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1770 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1771 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1772 cfmc      write(iout,*)'INTPRT=',intprt
1773 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1774 cfmc      write(iout,*)'IPRT=',iprt
1775 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1776 cfmc      write(iout,*)'IMAXTR=',imaxtr
1777 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1778 cfmc      write(iout,*)'MAXEVEN=',maxeven
1779 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1780 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1781 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1782 cfmc      write(iout,*)'INIMIN=',inimin
1783 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1784 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1785 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1786 cfmc      write(iout,*)'NTHREAD=',nthread
1787 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1788 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1789 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1790 cfmc      write(iout,*)'MAXPERT=',maxpert
1791 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1792 cfmc      write(iout,*)'IRMSD=',irmsd
1793 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1794 cfmc      write(iout,*)'DENEMIN=',denemin
1795 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1796 cfmc      write(iout,*)'RCUT1S=',rcut1s
1797 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1798 cfmc      write(iout,*)'RCUT1E=',rcut1e
1799 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1800 cfmc      write(iout,*)'RCUT2S=',rcut2s
1801 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1802 cfmc      write(iout,*)'RCUT2E=',rcut2e
1803 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1804 cfmc      write(iout,*)'DPERT1=',d_pert1
1805 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1806 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1807 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1808 cfmc      write(iout,*)'DPERT2=',d_pert2
1809 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1810 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1811 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1812 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1813 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1814 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1815 cfmc      d_pert1=deg2rad*d_pert1
1816 cfmc      d_pert1a=deg2rad*d_pert1a
1817 cfmc      d_pert2=deg2rad*d_pert2
1818 cfmc      d_pert2a=deg2rad*d_pert2a
1819 cfmc      d_pert2b=deg2rad*d_pert2b
1820 cfmc      d_pert2c=deg2rad*d_pert2c
1821 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1822 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1823 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1824 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1825 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1826 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1827 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1828 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1829 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1830 cfmc      write(iout,*)'RCUTINI=',rcutini
1831 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1832 cfmc      write(iout,*)'GRAT=',grat
1833 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1834 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1835 cfmc
1836 cfmc      return
1837 cfmc      end
1838 c----------------------------------------------------------------------------
1839       subroutine mcmread
1840       implicit real*8 (a-h,o-z)
1841       include 'DIMENSIONS'
1842       include 'COMMON.MCM'
1843       include 'COMMON.MCE'
1844       include 'COMMON.IOUNITS'
1845       character*80 ucase
1846       character*320 mcmcard
1847       call card_concat(mcmcard)
1848       call readi(mcmcard,'MAXACC',maxacc,100)
1849       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1850       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1851       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1852       call readi(mcmcard,'MAXREPM',maxrepm,200)
1853       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1854       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1855       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1856       call reada(mcmcard,'E_UP',e_up,5.0D0)
1857       call reada(mcmcard,'DELTE',delte,0.1D0)
1858       call readi(mcmcard,'NSWEEP',nsweep,5)
1859       call readi(mcmcard,'NSTEPH',nsteph,0)
1860       call readi(mcmcard,'NSTEPC',nstepc,0)
1861       call reada(mcmcard,'TMIN',tmin,298.0D0)
1862       call reada(mcmcard,'TMAX',tmax,298.0D0)
1863       call readi(mcmcard,'NWINDOW',nwindow,0)
1864       call readi(mcmcard,'PRINT_MC',print_mc,0)
1865       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1866       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1867       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1868       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1869       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1870       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1871       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1872       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1873       if (nwindow.gt.0) then
1874         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1875         do i=1,nwindow
1876           winlen(i)=winend(i)-winstart(i)+1
1877         enddo
1878       endif
1879       if (tmax.lt.tmin) tmax=tmin
1880       if (tmax.eq.tmin) then
1881         nstepc=0
1882         nsteph=0
1883       endif
1884       if (nstepc.gt.0 .and. nsteph.gt.0) then
1885         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1886         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1887       endif
1888 C Probabilities of different move types
1889       sumpro_type(0)=0.0D0
1890       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1891       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1892       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1893       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1894       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1895       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1896       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1897       do i=1,MaxMoveType
1898         print *,'i',i,' sumprotype',sumpro_type(i)
1899         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1900         print *,'i',i,' sumprotype',sumpro_type(i)
1901       enddo
1902       return
1903       end
1904 c----------------------------------------------------------------------------
1905       subroutine read_minim
1906       implicit real*8 (a-h,o-z)
1907       include 'DIMENSIONS'
1908       include 'COMMON.MINIM'
1909       include 'COMMON.IOUNITS'
1910       character*80 ucase
1911       character*320 minimcard
1912       call card_concat(minimcard)
1913       call readi(minimcard,'MAXMIN',maxmin,2000)
1914       call readi(minimcard,'MAXFUN',maxfun,5000)
1915       call readi(minimcard,'MINMIN',minmin,maxmin)
1916       call readi(minimcard,'MINFUN',minfun,maxmin)
1917       call reada(minimcard,'TOLF',tolf,1.0D-2)
1918       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1919       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1920       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1921       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1922       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1923      &         'Options in energy minimization:'
1924       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1925      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1926      & 'MinMin:',MinMin,' MinFun:',MinFun,
1927      & ' TolF:',TolF,' RTolF:',RTolF
1928       return
1929       end
1930 c----------------------------------------------------------------------------
1931       subroutine read_angles(kanal,*)
1932       implicit real*8 (a-h,o-z)
1933       include 'DIMENSIONS'
1934       include 'COMMON.GEO'
1935       include 'COMMON.VAR'
1936       include 'COMMON.CHAIN'
1937       include 'COMMON.IOUNITS'
1938       include 'COMMON.CONTROL'
1939 c Read angles from input
1940 c
1941        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1942        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1943        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1944        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1945
1946        do i=1,nres
1947 c 9/7/01 avoid 180 deg valence angle
1948         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1949 c
1950         theta(i)=deg2rad*theta(i)
1951         phi(i)=deg2rad*phi(i)
1952         alph(i)=deg2rad*alph(i)
1953         omeg(i)=deg2rad*omeg(i)
1954        enddo
1955       return
1956    10 return1
1957       end
1958 c----------------------------------------------------------------------------
1959       subroutine reada(rekord,lancuch,wartosc,default)
1960       implicit none
1961       character*(*) rekord,lancuch
1962       double precision wartosc,default
1963       integer ilen,iread
1964       external ilen
1965       iread=index(rekord,lancuch)
1966       if (iread.eq.0) then
1967         wartosc=default
1968         return
1969       endif
1970       iread=iread+ilen(lancuch)+1
1971       read (rekord(iread:),*,err=10,end=10) wartosc
1972       return
1973   10  wartosc=default
1974       return
1975       end
1976 c----------------------------------------------------------------------------
1977       subroutine readi(rekord,lancuch,wartosc,default)
1978       implicit none
1979       character*(*) rekord,lancuch
1980       integer wartosc,default
1981       integer ilen,iread
1982       external ilen
1983       iread=index(rekord,lancuch)
1984       if (iread.eq.0) then
1985         wartosc=default
1986         return
1987       endif
1988       iread=iread+ilen(lancuch)+1
1989       read (rekord(iread:),*,err=10,end=10) wartosc
1990       return
1991   10  wartosc=default
1992       return
1993       end
1994 c----------------------------------------------------------------------------
1995       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1996       implicit none
1997       integer dim,i
1998       integer tablica(dim),default
1999       character*(*) rekord,lancuch
2000       character*80 aux
2001       integer ilen,iread
2002       external ilen
2003       do i=1,dim
2004         tablica(i)=default
2005       enddo
2006       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2007       if (iread.eq.0) return
2008       iread=iread+ilen(lancuch)+1
2009       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2010    10 return
2011       end
2012 c----------------------------------------------------------------------------
2013       subroutine multreada(rekord,lancuch,tablica,dim,default)
2014       implicit none
2015       integer dim,i
2016       double precision tablica(dim),default
2017       character*(*) rekord,lancuch
2018       character*80 aux
2019       integer ilen,iread
2020       external ilen
2021       do i=1,dim
2022         tablica(i)=default
2023       enddo
2024       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2025       if (iread.eq.0) return
2026       iread=iread+ilen(lancuch)+1
2027       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2028    10 return
2029       end
2030 c----------------------------------------------------------------------------
2031       subroutine openunits
2032       implicit real*8 (a-h,o-z)
2033       include 'DIMENSIONS'
2034 #ifdef MPI
2035       include 'mpif.h'
2036       character*16 form,nodename
2037       integer nodelen
2038 #endif
2039       include 'COMMON.SETUP'
2040       include 'COMMON.IOUNITS'
2041       include 'COMMON.CONTROL'
2042       integer lenpre,lenpot,ilen,lentmp
2043       external ilen
2044       character*3 out1file_text,ucase
2045       character*3 ll
2046       external ucase
2047       tmpdir=""
2048       out1file_text="YES"
2049 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2050       call getenv_loc("PREFIX",prefix)
2051       pref_orig = prefix
2052       call getenv_loc("POT",pot)
2053       call getenv_loc("DIRTMP",tmpdir)
2054       call getenv_loc("CURDIR",curdir)
2055       call getenv_loc("OUT1FILE",out1file_text)
2056 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2057       out1file_text=ucase(out1file_text)
2058       if (out1file_text(1:1).eq."Y") then
2059         out1file=.true.
2060       else
2061         out1file=fg_rank.gt.0
2062       endif
2063       lenpre=ilen(prefix)
2064       lenpot=ilen(pot)
2065       lentmp=ilen(tmpdir)
2066       if (lentmp.gt.0) then
2067           write (*,'(80(1h!))')
2068           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2069           write (*,'(80(1h!))')
2070           write (*,*)"All output files will be on node /tmp directory."
2071 #ifdef MPI
2072         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2073         if (me.eq.king) then
2074           write (*,*) "The master node is ",nodename
2075         else if (fg_rank.eq.0) then
2076           write (*,*) "I am the CG slave node ",nodename
2077         else
2078           write (*,*) "I am the FG slave node ",nodename
2079         endif
2080 #endif
2081         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2082         lenpre = lentmp+lenpre+1
2083       endif
2084       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2085 C Get the names and open the input files
2086 #if defined(WINIFL) || defined(WINPGI)
2087       open(1,file=pref_orig(:ilen(pref_orig))//
2088      &  '.inp',status='old',readonly,shared)
2089        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2090 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2091 C Get parameter filenames and open the parameter files.
2092       call getenv_loc('BONDPAR',bondname)
2093       open (ibond,file=bondname,status='old',readonly,shared)
2094       call getenv_loc('THETPAR',thetname)
2095       open (ithep,file=thetname,status='old',readonly,shared)
2096       call getenv_loc('ROTPAR',rotname)
2097       open (irotam,file=rotname,status='old',readonly,shared)
2098       call getenv_loc('TORPAR',torname)
2099       open (itorp,file=torname,status='old',readonly,shared)
2100       call getenv_loc('TORDPAR',tordname)
2101       open (itordp,file=tordname,status='old',readonly,shared)
2102       call getenv_loc('FOURIER',fouriername)
2103       open (ifourier,file=fouriername,status='old',readonly,shared)
2104       call getenv_loc('ELEPAR',elename)
2105       open (ielep,file=elename,status='old',readonly,shared)
2106       call getenv_loc('SIDEPAR',sidename)
2107       open (isidep,file=sidename,status='old',readonly,shared)
2108       call getenv_loc('LIPTRANPAR',liptranname)
2109       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2110 #elif (defined CRAY) || (defined AIX)
2111       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2112      &  action='read')
2113 c      print *,"Processor",myrank," opened file 1"
2114       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2115 c      print *,"Processor",myrank," opened file 9"
2116 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2117 C Get parameter filenames and open the parameter files.
2118       call getenv_loc('BONDPAR',bondname)
2119       open (ibond,file=bondname,status='old',action='read')
2120 c      print *,"Processor",myrank," opened file IBOND"
2121       call getenv_loc('THETPAR',thetname)
2122       open (ithep,file=thetname,status='old',action='read')
2123 #ifndef CRYST_THETA
2124       call getenv_loc('THETPARPDB',thetname_pdb)
2125       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2126 #endif
2127 c      print *,"Processor",myrank," opened file ITHEP"
2128       call getenv_loc('ROTPAR',rotname)
2129       open (irotam,file=rotname,status='old',action='read')
2130 #ifndef CRYST_SC
2131       call getenv_loc('ROTPARPDB',rotname_pdb)
2132       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2133 #endif
2134 c      print *,"Processor",myrank," opened file IROTAM"
2135       call getenv_loc('TORPAR',torname)
2136       open (itorp,file=torname,status='old',action='read')
2137 c      print *,"Processor",myrank," opened file ITORP"
2138       call getenv_loc('TORDPAR',tordname)
2139       open (itordp,file=tordname,status='old',action='read')
2140 c      print *,"Processor",myrank," opened file ITORDP"
2141       call getenv_loc('SCCORPAR',sccorname)
2142       open (isccor,file=sccorname,status='old',action='read')
2143 c      print *,"Processor",myrank," opened file ISCCOR"
2144       call getenv_loc('FOURIER',fouriername)
2145       open (ifourier,file=fouriername,status='old',action='read')
2146 c      print *,"Processor",myrank," opened file IFOURIER"
2147       call getenv_loc('ELEPAR',elename)
2148       open (ielep,file=elename,status='old',action='read')
2149 c      print *,"Processor",myrank," opened file IELEP"
2150       call getenv_loc('SIDEPAR',sidename)
2151       open (isidep,file=sidename,status='old',action='read')
2152       call getenv_loc('LIPTRANPAR',liptranname)
2153       open (iliptranpar,file=liptranname,status='old',action='read')
2154 c      print *,"Processor",myrank," opened file ISIDEP"
2155 c      print *,"Processor",myrank," opened parameter files"
2156 #elif (defined G77)
2157       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2158       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2159 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2160 C Get parameter filenames and open the parameter files.
2161       call getenv_loc('BONDPAR',bondname)
2162       open (ibond,file=bondname,status='old')
2163       call getenv_loc('THETPAR',thetname)
2164       open (ithep,file=thetname,status='old')
2165 #ifndef CRYST_THETA
2166       call getenv_loc('THETPARPDB',thetname_pdb)
2167       open (ithep_pdb,file=thetname_pdb,status='old')
2168 #endif
2169       call getenv_loc('ROTPAR',rotname)
2170       open (irotam,file=rotname,status='old')
2171 #ifndef CRYST_SC
2172       call getenv_loc('ROTPARPDB',rotname_pdb)
2173       open (irotam_pdb,file=rotname_pdb,status='old')
2174 #endif
2175       call getenv_loc('TORPAR',torname)
2176       open (itorp,file=torname,status='old')
2177       call getenv_loc('TORDPAR',tordname)
2178       open (itordp,file=tordname,status='old')
2179       call getenv_loc('SCCORPAR',sccorname)
2180       open (isccor,file=sccorname,status='old')
2181       call getenv_loc('FOURIER',fouriername)
2182       open (ifourier,file=fouriername,status='old')
2183       call getenv_loc('ELEPAR',elename)
2184       open (ielep,file=elename,status='old')
2185       call getenv_loc('SIDEPAR',sidename)
2186       open (isidep,file=sidename,status='old')
2187       call getenv_loc('LIPTRANPAR',liptranname)
2188       open (iliptranpar,file=liptranname,status='old')
2189 #else
2190       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2191      &  readonly)
2192        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2193 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2194 C Get parameter filenames and open the parameter files.
2195       call getenv_loc('BONDPAR',bondname)
2196       open (ibond,file=bondname,status='old',readonly)
2197       call getenv_loc('THETPAR',thetname)
2198       open (ithep,file=thetname,status='old',readonly)
2199       call getenv_loc('ROTPAR',rotname)
2200       open (irotam,file=rotname,status='old',readonly)
2201       call getenv_loc('TORPAR',torname)
2202       open (itorp,file=torname,status='old',readonly)
2203       call getenv_loc('TORDPAR',tordname)
2204       open (itordp,file=tordname,status='old',readonly)
2205       call getenv_loc('SCCORPAR',sccorname)
2206       open (isccor,file=sccorname,status='old',readonly)
2207 #ifndef CRYST_THETA
2208       call getenv_loc('THETPARPDB',thetname_pdb)
2209       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2210 #endif
2211       call getenv_loc('FOURIER',fouriername)
2212       open (ifourier,file=fouriername,status='old',readonly)
2213       call getenv_loc('ELEPAR',elename)
2214       open (ielep,file=elename,status='old',readonly)
2215       call getenv_loc('SIDEPAR',sidename)
2216       open (isidep,file=sidename,status='old',readonly)
2217       call getenv_loc('LIPTRANPAR',liptranname)
2218       open (iliptranpar,file=liptranname,status='old',action='read')
2219 #ifndef CRYST_SC
2220       call getenv_loc('ROTPARPDB',rotname_pdb)
2221       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2222 #endif
2223 #endif
2224 #ifdef TUBE
2225       call getenv_loc('TUBEPAR',tubename)
2226 #if defined(WINIFL) || defined(WINPGI)
2227       open (itube,file=tubename,status='old',readonly,shared)
2228 #elif (defined CRAY)  || (defined AIX)
2229       open (itube,file=tubename,status='old',action='read')
2230 #elif (defined G77)
2231       open (itube,file=tubename,status='old')
2232 #else
2233       open (itube,file=tubename,status='old',readonly)
2234 #endif
2235 #endif
2236 #ifndef OLDSCP
2237 C
2238 C 8/9/01 In the newest version SCp interaction constants are read from a file
2239 C Use -DOLDSCP to use hard-coded constants instead.
2240 C
2241       call getenv_loc('SCPPAR',scpname)
2242 #if defined(WINIFL) || defined(WINPGI)
2243       open (iscpp,file=scpname,status='old',readonly,shared)
2244 #elif (defined CRAY)  || (defined AIX)
2245       open (iscpp,file=scpname,status='old',action='read')
2246 #elif (defined G77)
2247       open (iscpp,file=scpname,status='old')
2248 #else
2249       open (iscpp,file=scpname,status='old',readonly)
2250 #endif
2251 #endif
2252       call getenv_loc('PATTERN',patname)
2253 #if defined(WINIFL) || defined(WINPGI)
2254       open (icbase,file=patname,status='old',readonly,shared)
2255 #elif (defined CRAY)  || (defined AIX)
2256       open (icbase,file=patname,status='old',action='read')
2257 #elif (defined G77)
2258       open (icbase,file=patname,status='old')
2259 #else
2260       open (icbase,file=patname,status='old',readonly)
2261 #endif
2262 #ifdef MPI
2263 C Open output file only for CG processes
2264 c      print *,"Processor",myrank," fg_rank",fg_rank
2265       if (fg_rank.eq.0) then
2266
2267       if (nodes.eq.1) then
2268         npos=3
2269       else
2270         npos = dlog10(dfloat(nodes-1))+1
2271       endif
2272       if (npos.lt.3) npos=3
2273       write (liczba,'(i1)') npos
2274       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2275      &  //')'
2276       write (liczba,form) me
2277       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2278      &  liczba(:ilen(liczba))
2279       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2280      &  //'.int'
2281       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2282      &  //'.pdb'
2283       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2284      &  liczba(:ilen(liczba))//'.mol2'
2285       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2286      &  liczba(:ilen(liczba))//'.stat'
2287       if (lentmp.gt.0)
2288      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2289      &      //liczba(:ilen(liczba))//'.stat')
2290       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2291      &  //'.rst'
2292       if(usampl) then
2293           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2294      & liczba(:ilen(liczba))//'.const'
2295       endif
2296
2297       endif
2298 #else
2299       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2300       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2301       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2302       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2303       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2304       if (lentmp.gt.0)
2305      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2306      &    //'.stat')
2307       rest2name=prefix(:ilen(prefix))//'.rst'
2308       if(usampl) then
2309          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2310       endif
2311 #endif
2312 #if defined(AIX) || defined(PGI) || defined(CRAY)
2313       if (me.eq.king .or. .not. out1file) then
2314          open(iout,file=outname,status='unknown')
2315       else
2316          open(iout,file="/dev/null",status="unknown")
2317       endif
2318 c#define DEBUG
2319 #ifdef DEBUG
2320       if (fg_rank.gt.0) then
2321         write (liczba,'(i3.3)') myrank/nfgtasks
2322         write (ll,'(bz,i3.3)') fg_rank
2323         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2324      &   status='unknown')
2325       endif
2326 #endif
2327 c#undef DEBUG
2328       if(me.eq.king) then
2329        open(igeom,file=intname,status='unknown',position='append')
2330        open(ipdb,file=pdbname,status='unknown')
2331        open(imol2,file=mol2name,status='unknown')
2332        open(istat,file=statname,status='unknown',position='append')
2333       else
2334 c1out       open(iout,file=outname,status='unknown')
2335       endif
2336 #else
2337       if (me.eq.king .or. .not.out1file)
2338      &    open(iout,file=outname,status='unknown')
2339 #ifdef DEBUG
2340       if (fg_rank.gt.0) then
2341         write (liczba,'(i3.3)') myrank/nfgtasks
2342         write (ll,'(bz,i3.3)') fg_rank
2343         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2344      &   status='unknown')
2345       endif
2346 #endif
2347       if(me.eq.king) then
2348        open(igeom,file=intname,status='unknown',access='append')
2349        open(ipdb,file=pdbname,status='unknown')
2350        open(imol2,file=mol2name,status='unknown')
2351        open(istat,file=statname,status='unknown',access='append')
2352       else
2353 c1out       open(iout,file=outname,status='unknown')
2354       endif
2355 #endif
2356       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2357       csa_seed=prefix(:lenpre)//'.CSA.seed'
2358       csa_history=prefix(:lenpre)//'.CSA.history'
2359       csa_bank=prefix(:lenpre)//'.CSA.bank'
2360       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2361       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2362       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2363 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2364       csa_int=prefix(:lenpre)//'.int'
2365       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2366       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2367       csa_in=prefix(:lenpre)//'.CSA.in'
2368 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2369 C Write file names
2370       if (me.eq.king)then
2371       write (iout,'(80(1h-))')
2372       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2373       write (iout,'(80(1h-))')
2374       write (iout,*) "Input file                      : ",
2375      &  pref_orig(:ilen(pref_orig))//'.inp'
2376       write (iout,*) "Output file                     : ",
2377      &  outname(:ilen(outname))
2378       write (iout,*)
2379       write (iout,*) "Sidechain potential file        : ",
2380      &  sidename(:ilen(sidename))
2381 #ifndef OLDSCP
2382       write (iout,*) "SCp potential file              : ",
2383      &  scpname(:ilen(scpname))
2384 #endif
2385       write (iout,*) "Electrostatic potential file    : ",
2386      &  elename(:ilen(elename))
2387       write (iout,*) "Cumulant coefficient file       : ",
2388      &  fouriername(:ilen(fouriername))
2389       write (iout,*) "Torsional parameter file        : ",
2390      &  torname(:ilen(torname))
2391       write (iout,*) "Double torsional parameter file : ",
2392      &  tordname(:ilen(tordname))
2393       write (iout,*) "SCCOR parameter file : ",
2394      &  sccorname(:ilen(sccorname))
2395       write (iout,*) "Bond & inertia constant file    : ",
2396      &  bondname(:ilen(bondname))
2397       write (iout,*) "Bending-torsion parameter file  : ",
2398      &  thetname(:ilen(thetname))
2399       write (iout,*) "Rotamer parameter file          : ",
2400      &  rotname(:ilen(rotname))
2401       write (iout,*) "Thetpdb parameter file          : ",
2402      &  thetname_pdb(:ilen(thetname_pdb))
2403       write (iout,*) "Threading database              : ",
2404      &  patname(:ilen(patname))
2405       if (lentmp.ne.0)
2406      &write (iout,*)" DIRTMP                          : ",
2407      &  tmpdir(:lentmp)
2408       write (iout,'(80(1h-))')
2409       endif
2410       return
2411       end
2412 c----------------------------------------------------------------------------
2413       subroutine card_concat(card)
2414       implicit real*8 (a-h,o-z)
2415       include 'DIMENSIONS'
2416       include 'COMMON.IOUNITS'
2417       character*(*) card
2418       character*80 karta,ucase
2419       external ilen
2420       read (inp,'(a)') karta
2421       karta=ucase(karta)
2422       card=' '
2423       do while (karta(80:80).eq.'&')
2424         card=card(:ilen(card)+1)//karta(:79)
2425         read (inp,'(a)') karta
2426         karta=ucase(karta)
2427       enddo
2428       card=card(:ilen(card)+1)//karta
2429       return
2430       end
2431 c----------------------------------------------------------------------------------
2432       subroutine readrst
2433       implicit real*8 (a-h,o-z)
2434       include 'DIMENSIONS'
2435       include 'COMMON.CHAIN'
2436       include 'COMMON.IOUNITS'
2437       include 'COMMON.CONTROL'
2438       include 'COMMON.LAGRANGE'
2439       open(irest2,file=rest2name,status='unknown')
2440       read(irest2,*) totT,EK,potE,totE,t_bath
2441       totTafm=totT
2442       do i=1,2*nres
2443          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2444       enddo
2445       do i=1,2*nres
2446          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2447       enddo
2448       if(usampl) then
2449              read (irest2,*) iset
2450       endif
2451       close(irest2)
2452       return
2453       end
2454 c---------------------------------------------------------------------------------
2455       subroutine read_fragments
2456       implicit real*8 (a-h,o-z)
2457       include 'DIMENSIONS'
2458 #ifdef MPI
2459       include 'mpif.h'
2460 #endif
2461       include 'COMMON.SETUP'
2462       include 'COMMON.CHAIN'
2463       include 'COMMON.IOUNITS'
2464       include 'COMMON.QRESTR'
2465       include 'COMMON.CONTROL'
2466       read(inp,*) nset,nfrag,npair,nfrag_back
2467       loc_qlike=(nfrag_back.lt.0)
2468       nfrag_back=iabs(nfrag_back)
2469       if(me.eq.king.or..not.out1file)
2470      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2471      &  " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2472       do iset=1,nset
2473          read(inp,*) mset(iset)
2474        do i=1,nfrag
2475          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2476      &     qinfrag(i,iset)
2477          if(me.eq.king.or..not.out1file)
2478      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2479      &     ifrag(2,i,iset), qinfrag(i,iset)
2480        enddo
2481        do i=1,npair
2482         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2483      &    qinpair(i,iset)
2484         if(me.eq.king.or..not.out1file)
2485      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2486      &    ipair(2,i,iset), qinpair(i,iset)
2487        enddo
2488        if (loc_qlike) then
2489        do i=1,nfrag_back
2490         read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2491      &     wfrag_back(2,i,iset),qin_back(2,i,iset),
2492      &     wfrag_back(3,i,iset),qin_back(3,i,iset),
2493      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2494         if(me.eq.king.or..not.out1file)
2495      &   write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2496      &   wfrag_back(2,i,iset),qin_back(3,i,iset),
2497      &   wfrag_back(3,i,iset),qin_back(3,i,iset),
2498      &   ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2499        enddo
2500        else
2501        do i=1,nfrag_back
2502         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2503      &     wfrag_back(3,i,iset),
2504      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2505         if(me.eq.king.or..not.out1file)
2506      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2507      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2508        enddo
2509        endif
2510       enddo
2511       return
2512       end
2513 C---------------------------------------------------------------------------
2514       subroutine read_afminp
2515             implicit real*8 (a-h,o-z)
2516       include 'DIMENSIONS'
2517 #ifdef MPI
2518       include 'mpif.h'
2519 #endif
2520       include 'COMMON.SETUP'
2521       include 'COMMON.CONTROL'
2522       include 'COMMON.CHAIN'
2523       include 'COMMON.IOUNITS'
2524       include 'COMMON.SBRIDGE'
2525       character*320 afmcard
2526       print *, "wchodze"
2527       call card_concat(afmcard)
2528       call readi(afmcard,"BEG",afmbeg,0)
2529       call readi(afmcard,"END",afmend,0)
2530       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2531       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2532       print *,'FORCE=' ,forceAFMconst
2533 CCCC NOW PROPERTIES FOR AFM
2534        distafminit=0.0d0
2535        do i=1,3
2536         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2537        enddo
2538         distafminit=dsqrt(distafminit)
2539         print *,'initdist',distafminit
2540       return
2541       end
2542 c-------------------------------------------------------------------------------
2543       subroutine read_saxs_constr
2544       implicit real*8 (a-h,o-z)
2545       include 'DIMENSIONS'
2546 #ifdef MPI
2547       include 'mpif.h'
2548 #endif
2549       include 'COMMON.SETUP'
2550       include 'COMMON.CONTROL'
2551       include 'COMMON.CHAIN'
2552       include 'COMMON.IOUNITS'
2553       include 'COMMON.SBRIDGE'
2554       include 'COMMON.SAXS'
2555       double precision cm(3)
2556 c      read(inp,*) nsaxs
2557       write (iout,*) "Calling read_saxs nsaxs",nsaxs
2558       call flush(iout)
2559       if (saxs_mode.eq.0) then
2560 c SAXS distance distribution
2561       do i=1,nsaxs
2562         read(inp,*) distsaxs(i),Psaxs(i)
2563       enddo
2564       Cnorm = 0.0d0
2565       do i=1,nsaxs
2566         Cnorm = Cnorm + Psaxs(i)
2567       enddo
2568       write (iout,*) "Cnorm",Cnorm
2569       do i=1,nsaxs
2570         Psaxs(i)=Psaxs(i)/Cnorm
2571       enddo
2572       write (iout,*) "Normalized distance distribution from SAXS"
2573       do i=1,nsaxs
2574         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2575       enddo
2576       Wsaxs0=0.0d0
2577       do i=1,nsaxs
2578         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2579       enddo
2580       write (iout,*) "Wsaxs0",Wsaxs0
2581       else
2582 c SAXS "spheres".
2583       do i=1,nsaxs
2584         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2585       enddo
2586       do j=1,3
2587         cm(j)=0.0d0
2588       enddo
2589       do i=1,nsaxs
2590         do j=1,3
2591           cm(j)=cm(j)+Csaxs(j,i)
2592         enddo
2593       enddo
2594       do j=1,3
2595         cm(j)=cm(j)/nsaxs
2596       enddo
2597       do i=1,nsaxs
2598         do j=1,3
2599           Csaxs(j,i)=Csaxs(j,i)-cm(j)
2600         enddo
2601       enddo
2602       write (iout,*) "SAXS sphere coordinates"
2603       do i=1,nsaxs
2604         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2605       enddo
2606       endif
2607       return
2608       end
2609
2610 c-------------------------------------------------------------------------------
2611       subroutine read_dist_constr
2612       implicit real*8 (a-h,o-z)
2613       include 'DIMENSIONS'
2614 #ifdef MPI
2615       include 'mpif.h'
2616 #endif
2617       include 'COMMON.SETUP'
2618       include 'COMMON.CONTROL'
2619       include 'COMMON.CHAIN'
2620       include 'COMMON.IOUNITS'
2621       include 'COMMON.SBRIDGE'
2622       include 'COMMON.INTERACT'
2623       integer ifrag_(2,100),ipair_(2,1000)
2624       double precision wfrag_(100),wpair_(1000)
2625       character*5000 controlcard
2626       logical normalize,next
2627       integer restr_type
2628       double precision scal_bfac
2629       double precision xlink(4,0:4) /
2630 c           a          b       c     sigma
2631      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
2632      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
2633      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
2634      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
2635      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
2636 c      print *, "WCHODZE"
2637       write (iout,*) "Calling read_dist_constr"
2638 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2639 c      call flush(iout)
2640       restr_on_coord=.false.
2641       next=.true.
2642
2643       npeak=0
2644       ipeak=0
2645       nhpb_peak=0
2646
2647       DO WHILE (next)
2648
2649       call card_concat(controlcard)
2650       next = index(controlcard,"NEXT").gt.0
2651       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2652       write (iout,*) "restr_type",restr_type
2653       call readi(controlcard,"NFRAG",nfrag_,0)
2654       call readi(controlcard,"NFRAG",nfrag_,0)
2655       call readi(controlcard,"NPAIR",npair_,0)
2656       call readi(controlcard,"NDIST",ndist_,0)
2657       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2658       call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2659       if (restr_type.eq.10)
2660      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2661       if (restr_type.eq.12)
2662      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2663       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2664       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2665       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2666       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2667       normalize = index(controlcard,"NORMALIZE").gt.0
2668       write (iout,*) "WBOLTZD",wboltzd
2669       write (iout,*) "SCAL_PEAK",scal_peak
2670       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2671       write (iout,*) "IFRAG"
2672       do i=1,nfrag_
2673         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2674       enddo
2675       write (iout,*) "IPAIR"
2676       do i=1,npair_
2677         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2678       enddo
2679       if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2680      & write (iout,*)
2681      &   "Distance restraints as generated from reference structure"
2682       do i=1,nfrag_
2683         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2684         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2685      &    ifrag_(2,i)=nstart_sup+nsup-1
2686 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2687 c        call flush(iout)
2688         if (wfrag_(i).eq.0.0d0) cycle
2689         do j=ifrag_(1,i),ifrag_(2,i)-1
2690           do k=j+1,ifrag_(2,i)
2691 c            write (iout,*) "j",j," k",k
2692             ddjk=dist(j,k)
2693             if (restr_type.eq.1) then
2694               nhpb=nhpb+1
2695               irestr_type(nhpb)=1
2696               ihpb(nhpb)=j
2697               jhpb(nhpb)=k
2698               dhpb(nhpb)=ddjk
2699               forcon(nhpb)=wfrag_(i)
2700             else if (constr_dist.eq.2) then
2701               if (ddjk.le.dist_cut) then
2702                 nhpb=nhpb+1
2703                 irestr_type(nhpb)=1
2704                 ihpb(nhpb)=j
2705                 jhpb(nhpb)=k
2706                 dhpb(nhpb)=ddjk
2707                 forcon(nhpb)=wfrag_(i)
2708               endif
2709             else if (restr_type.eq.3) then
2710               nhpb=nhpb+1
2711               irestr_type(nhpb)=1
2712               ihpb(nhpb)=j
2713               jhpb(nhpb)=k
2714               dhpb(nhpb)=ddjk
2715               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2716             endif
2717 #ifdef MPI
2718             if (.not.out1file .or. me.eq.king)
2719      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2720      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2721 #else
2722             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2723      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2724 #endif
2725           enddo
2726         enddo
2727       enddo
2728       do i=1,npair_
2729         if (wpair_(i).eq.0.0d0) cycle
2730         ii = ipair_(1,i)
2731         jj = ipair_(2,i)
2732         if (ii.gt.jj) then
2733           itemp=ii
2734           ii=jj
2735           jj=itemp
2736         endif
2737         do j=ifrag_(1,ii),ifrag_(2,ii)
2738           do k=ifrag_(1,jj),ifrag_(2,jj)
2739             ddjk=dist(j,k)
2740             if (restr_type.eq.1) then
2741               nhpb=nhpb+1
2742               irestr_type(nhpb)=1
2743               ihpb(nhpb)=j
2744               jhpb(nhpb)=k
2745               dhpb(nhpb)=ddjk
2746               forcon(nhpb)=wpair_(i)
2747             else if (constr_dist.eq.2) then
2748               if (ddjk.le.dist_cut) then
2749                 nhpb=nhpb+1
2750                 irestr_type(nhpb)=1
2751                 ihpb(nhpb)=j
2752                 jhpb(nhpb)=k
2753                 dhpb(nhpb)=ddjk
2754                 forcon(nhpb)=wpair_(i)
2755               endif
2756             else if (restr_type.eq.3) then
2757               nhpb=nhpb+1
2758               irestr_type(nhpb)=1
2759               ihpb(nhpb)=j
2760               jhpb(nhpb)=k
2761               dhpb(nhpb)=ddjk
2762               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2763             endif
2764 #ifdef MPI
2765             if (.not.out1file .or. me.eq.king)
2766      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2767      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2768 #else
2769             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2770      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2771 #endif
2772           enddo
2773         enddo
2774       enddo
2775
2776 c      print *,ndist_
2777       write (iout,*) "Distance restraints as read from input"
2778       do i=1,ndist_
2779         if (restr_type.eq.12) then
2780           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2781      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2782      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2783      &    fordepth_peak(nhpb_peak+1),npeak
2784 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2785 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2786 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2787 c     &    fordepth_peak(nhpb_peak+1),npeak
2788           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2789      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2790           nhpb_peak=nhpb_peak+1
2791           irestr_type_peak(nhpb_peak)=12
2792           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2793           ipeak(2,npeak)=i
2794 #ifdef MPI
2795           if (.not.out1file .or. me.eq.king)
2796      &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2797      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2798      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2799      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2800      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2801 #else
2802           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2803      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2804      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2805      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2806      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2807 #endif
2808           if (ibecarb_peak(nhpb_peak).eq.3) then
2809             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2810           else if (ibecarb_peak(nhpb_peak).eq.2) then
2811             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2812           else if (ibecarb_peak(nhpb_peak).eq.1) then
2813             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2814             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2815           endif
2816         else if (restr_type.eq.11) then
2817           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2818      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2819 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2820           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2821           nhpb=nhpb+1
2822           irestr_type(nhpb)=11
2823 #ifdef MPI
2824           if (.not.out1file .or. me.eq.king)
2825      &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2826      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2827      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2828 #else
2829           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2830      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2831      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2832 #endif
2833 c          if (ibecarb(nhpb).gt.0) then
2834 c            ihpb(nhpb)=ihpb(nhpb)+nres
2835 c            jhpb(nhpb)=jhpb(nhpb)+nres
2836 c          endif
2837          if (ibecarb(nhpb).eq.3) then
2838             ihpb(nhpb)=ihpb(nhpb)+nres
2839           else if (ibecarb(nhpb).eq.2) then
2840             ihpb(nhpb)=ihpb(nhpb)+nres
2841           else if (ibecarb(nhpb).eq.1) then
2842             ihpb(nhpb)=ihpb(nhpb)+nres
2843             jhpb(nhpb)=jhpb(nhpb)+nres
2844           endif
2845         else if (restr_type.eq.10) then
2846 c Cross-lonk Markov-like potential
2847           call card_concat(controlcard)
2848           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2849           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2850           ibecarb(nhpb+1)=0
2851           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2852           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2853           if (index(controlcard,"ZL").gt.0) then
2854             link_type=1
2855           else if (index(controlcard,"ADH").gt.0) then
2856             link_type=2
2857           else if (index(controlcard,"PDH").gt.0) then
2858             link_type=3
2859           else if (index(controlcard,"DSS").gt.0) then
2860             link_type=4
2861           else
2862             link_type=0
2863           endif
2864           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2865      &       xlink(1,link_type))
2866           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2867      &       xlink(2,link_type))
2868           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2869      &       xlink(3,link_type))
2870           call reada(controlcard,"SIGMA",forcon(nhpb+1),
2871      &       xlink(4,link_type))
2872           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2873 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2874 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2875           if (forcon(nhpb+1).le.0.0d0 .or.
2876      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2877           nhpb=nhpb+1
2878           irestr_type(nhpb)=10
2879           if (ibecarb(nhpb).eq.3) then
2880             jhpb(nhpb)=jhpb(nhpb)+nres
2881           else if (ibecarb(nhpb).eq.2) then
2882             ihpb(nhpb)=ihpb(nhpb)+nres
2883           else if (ibecarb(nhpb).eq.1) then
2884             ihpb(nhpb)=ihpb(nhpb)+nres
2885             jhpb(nhpb)=jhpb(nhpb)+nres
2886           endif
2887 #ifdef MPI
2888           if (.not.out1file .or. me.eq.king)
2889      &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2890      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2891      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2892      &     irestr_type(nhpb)
2893 #else
2894           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2895      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2896      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2897      &     irestr_type(nhpb)
2898 #endif
2899         else
2900 C        print *,"in else"
2901           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2902      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2903           if (forcon(nhpb+1).gt.0.0d0) then
2904           nhpb=nhpb+1
2905           if (dhpb1(nhpb).eq.0.0d0) then
2906             irestr_type(nhpb)=1
2907           else
2908             irestr_type(nhpb)=2
2909           endif
2910           if (ibecarb(nhpb).eq.3) then
2911             jhpb(nhpb)=jhpb(nhpb)+nres
2912           else if (ibecarb(nhpb).eq.2) then
2913             ihpb(nhpb)=ihpb(nhpb)+nres
2914           else if (ibecarb(nhpb).eq.1) then
2915             ihpb(nhpb)=ihpb(nhpb)+nres
2916             jhpb(nhpb)=jhpb(nhpb)+nres
2917           endif
2918           if (dhpb(nhpb).eq.0.0d0)
2919      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2920         endif
2921 #ifdef MPI
2922           if (.not.out1file .or. me.eq.king)
2923      &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2924      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2925 #else
2926           write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2927      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2928 #endif
2929         endif
2930 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2931 C        if (forcon(nhpb+1).gt.0.0d0) then
2932 C          nhpb=nhpb+1
2933 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2934       enddo
2935
2936       if (restr_type.eq.4) then
2937         write (iout,*) "The BFAC array"
2938         do i=nnt,nct
2939           write (iout,'(i5,f10.5)') i,bfac(i)
2940         enddo
2941         do i=nnt,nct
2942           if (itype(i).eq.ntyp1) cycle
2943           do j=nnt,i-1
2944             if (itype(j).eq.ntyp1) cycle
2945             if (itype(i).eq.10) then
2946               iiend=0
2947             else
2948               iiend=1
2949             endif
2950             if (itype(j).eq.10) then
2951               jjend=0
2952             else
2953               jjend=1
2954             endif
2955             kk=0
2956             do ii=0,iiend
2957             do jj=0,jjend
2958             nhpb=nhpb+1
2959             irestr_type(nhpb)=1
2960             forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2961             irestr_type(nhpb)=1
2962             ibecarb(nhpb)=kk
2963             if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2964             ihpb(nhpb)=i+nres*ii
2965             jhpb(nhpb)=j+nres*jj
2966             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2967 #ifdef MPI
2968             if (.not.out1file .or. me.eq.king) then
2969             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2970      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2971      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2972      &       irestr_type(nhpb)
2973             endif
2974 #else
2975             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2976      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2977      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2978      &       irestr_type(nhpb)
2979 #endif
2980             kk=kk+1
2981           enddo
2982           enddo
2983           enddo
2984         enddo
2985       endif
2986
2987       if (restr_type.eq.5) then
2988         restr_on_coord=.true.
2989         do i=nnt,nct
2990           if (itype(i).eq.ntyp1) cycle
2991           bfac(i)=(scal_bfac/bfac(i))**2
2992         enddo
2993       endif
2994
2995       ENDDO ! next
2996
2997       fordepthmax=0.0d0
2998       if (normalize) then
2999         do i=nss+1,nhpb
3000           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3001      &      fordepthmax=fordepth(i)
3002         enddo
3003         do i=nss+1,nhpb
3004           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3005         enddo
3006       endif
3007       return
3008       end
3009 c-------------------------------------------------------------------------------
3010       subroutine read_constr_homology
3011
3012       include 'DIMENSIONS'
3013 #ifdef MPI
3014       include 'mpif.h'
3015 #endif
3016       include 'COMMON.SETUP'
3017       include 'COMMON.CONTROL'
3018       include 'COMMON.CHAIN'
3019       include 'COMMON.IOUNITS'
3020       include 'COMMON.HOMOLOGY'
3021       include 'COMMON.MD'
3022       include 'COMMON.GEO'
3023       include 'COMMON.INTERACT'
3024       include 'COMMON.NAMES'
3025 c
3026 c For new homol impl
3027 c
3028       include 'COMMON.VAR'
3029 c
3030
3031 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3032 c    &                 dist_cut
3033 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3034 c    &    sigma_odl_temp(maxres,maxres,max_template)
3035       character*2 kic2
3036       character*24 model_ki_dist, model_ki_angle
3037       character*500 controlcard
3038       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3039       integer ilen
3040       external ilen
3041       logical liiflag
3042 c
3043 c     FP - Nov. 2014 Temporary specifications for new vars
3044 c
3045       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3046      &                 rescore3_tmp
3047       double precision, dimension (max_template,maxres) :: rescore
3048       double precision, dimension (max_template,maxres) :: rescore2
3049       double precision, dimension (max_template,maxres) :: rescore3
3050       character*24 pdbfile,tpl_k_rescore
3051 c -----------------------------------------------------------------
3052 c Reading multiple PDB ref structures and calculation of retraints
3053 c not using pre-computed ones stored in files model_ki_{dist,angle}
3054 c FP (Nov., 2014)
3055 c -----------------------------------------------------------------
3056 c
3057 c
3058 c Alternative: reading from input
3059       call card_concat(controlcard)
3060       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3061       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3062       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3063       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3064       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3065       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3066       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3067       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3068       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3069       if(.not.read2sigma.and.start_from_model) then
3070           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3071      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3072           start_from_model=.false.
3073       endif
3074       if(start_from_model .and. (me.eq.king .or. .not. out1file))
3075      &    write(iout,*) 'START_FROM_MODELS is ON'
3076       if(start_from_model .and. rest) then
3077         if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3078          write(iout,*) 'START_FROM_MODELS is OFF'
3079          write(iout,*) 'remove restart keyword from input'
3080         endif
3081       endif
3082       if (homol_nset.gt.1)then
3083          call card_concat(controlcard)
3084          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3085          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3086           write(iout,*) "iset homology_weight "
3087           do i=1,homol_nset
3088            write(iout,*) i,waga_homology(i)
3089           enddo
3090          endif
3091          iset=mod(kolor,homol_nset)+1
3092       else
3093        iset=1
3094        waga_homology(1)=1.0
3095       endif
3096
3097 cd      write (iout,*) "nnt",nnt," nct",nct
3098 cd      call flush(iout)
3099
3100
3101       lim_odl=0
3102       lim_dih=0
3103 c
3104 c      write(iout,*) 'nnt=',nnt,'nct=',nct
3105 c
3106       do i = nnt,nct
3107         do k=1,constr_homology
3108          idomain(k,i)=0
3109         enddo
3110       enddo
3111
3112       ii=0
3113       do i = nnt,nct-2
3114         do j=i+2,nct
3115         ii=ii+1
3116         ii_in_use(ii)=0
3117         enddo
3118       enddo
3119
3120       if (read_homol_frag) then
3121         call read_klapaucjusz
3122       else
3123
3124       do k=1,constr_homology
3125
3126         read(inp,'(a)') pdbfile
3127         if(me.eq.king .or. .not. out1file)
3128      &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3129      &  pdbfile(:ilen(pdbfile))
3130         open(ipdbin,file=pdbfile,status='old',err=33)
3131         goto 34
3132   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3133      &  pdbfile(:ilen(pdbfile))
3134         stop
3135   34    continue
3136 c        print *,'Begin reading pdb data'
3137 c
3138 c Files containing res sim or local scores (former containing sigmas)
3139 c
3140
3141         write(kic2,'(bz,i2.2)') k
3142
3143         tpl_k_rescore="template"//kic2//".sco"
3144
3145         unres_pdb=.false.
3146         if (read2sigma) then
3147           call readpdb_template(k)
3148         else
3149           call readpdb
3150         endif
3151 c
3152 c     Distance restraints
3153 c
3154 c          ... --> odl(k,ii)
3155 C Copy the coordinates from reference coordinates (?)
3156         do i=1,2*nres
3157           do j=1,3
3158             c(j,i)=cref(j,i)
3159 c           write (iout,*) "c(",j,i,") =",c(j,i)
3160           enddo
3161         enddo
3162 c
3163 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3164 c
3165 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3166           open (ientin,file=tpl_k_rescore,status='old')
3167           if (nnt.gt.1) rescore(k,1)=0.0d0
3168           do irec=nnt,nct ! loop for reading res sim
3169             if (read2sigma) then
3170              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3171      &                                rescore3_tmp,idomain_tmp
3172              i_tmp=i_tmp+nnt-1
3173              idomain(k,i_tmp)=idomain_tmp
3174              rescore(k,i_tmp)=rescore_tmp
3175              rescore2(k,i_tmp)=rescore2_tmp
3176              rescore3(k,i_tmp)=rescore3_tmp
3177              if (.not. out1file .or. me.eq.king)
3178      &        write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3179      &                      i_tmp,rescore2_tmp,rescore_tmp,
3180      &                                rescore3_tmp,idomain_tmp
3181             else
3182              idomain(k,irec)=1
3183              read (ientin,*,end=1401) rescore_tmp
3184
3185 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3186              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3187 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3188             endif
3189           enddo
3190  1401   continue
3191         close (ientin)
3192         if (waga_dist.ne.0.0d0) then
3193           ii=0
3194           do i = nnt,nct-2
3195             do j=i+2,nct
3196
3197               x12=c(1,i)-c(1,j)
3198               y12=c(2,i)-c(2,j)
3199               z12=c(3,i)-c(3,j)
3200               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3201 c              write (iout,*) k,i,j,distal,dist2_cut
3202
3203             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3204      &            .and. distal.le.dist2_cut ) then
3205
3206               ii=ii+1
3207               ii_in_use(ii)=1
3208               l_homo(k,ii)=.true.
3209
3210 c             write (iout,*) "k",k
3211 c             write (iout,*) "i",i," j",j," constr_homology",
3212 c    &                       constr_homology
3213               ires_homo(ii)=i
3214               jres_homo(ii)=j
3215               odl(k,ii)=distal
3216               if (read2sigma) then
3217                 sigma_odl(k,ii)=0
3218                 do ik=i,j
3219                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3220                 enddo
3221                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3222                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3223      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3224               else
3225                 if (odl(k,ii).le.dist_cut) then
3226                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3227                 else
3228 #ifdef OLDSIGMA
3229                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3230      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3231 #else
3232                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3233      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3234 #endif
3235                 endif
3236               endif
3237               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3238             else
3239               ii=ii+1
3240               l_homo(k,ii)=.false.
3241             endif
3242             enddo
3243           enddo
3244         lim_odl=ii
3245         endif
3246 c
3247 c     Theta, dihedral and SC retraints
3248 c
3249         if (waga_angle.gt.0.0d0) then
3250 c         open (ientin,file=tpl_k_sigma_dih,status='old')
3251 c         do irec=1,maxres-3 ! loop for reading sigma_dih
3252 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3253 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3254 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3255 c    &                            sigma_dih(k,i+nnt-1)
3256 c         enddo
3257 c1402   continue
3258 c         close (ientin)
3259           do i = nnt+3,nct
3260             if (idomain(k,i).eq.0) then
3261                sigma_dih(k,i)=0.0
3262                cycle
3263             endif
3264             dih(k,i)=phiref(i) ! right?
3265 c           read (ientin,*) sigma_dih(k,i) ! original variant
3266 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
3267 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3268 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3269 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
3270 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
3271
3272             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3273      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3274 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3275 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
3276 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3277 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3278 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
3279 c   Instead of res sim other local measure of b/b str reliability possible
3280             if (sigma_dih(k,i).ne.0)
3281      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3282 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3283           enddo
3284           lim_dih=nct-nnt-2
3285         endif
3286
3287         if (waga_theta.gt.0.0d0) then
3288 c         open (ientin,file=tpl_k_sigma_theta,status='old')
3289 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3290 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3291 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3292 c    &                              sigma_theta(k,i+nnt-1)
3293 c         enddo
3294 c1403   continue
3295 c         close (ientin)
3296
3297           do i = nnt+2,nct ! right? without parallel.
3298 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
3299 c         do i=ithet_start,ithet_end ! with FG parallel.
3300              if (idomain(k,i).eq.0) then
3301               sigma_theta(k,i)=0.0
3302               cycle
3303              endif
3304              thetatpl(k,i)=thetaref(i)
3305 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3306 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
3307 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3308 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
3309 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
3310              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3311      &                        rescore(k,i-2))/3.0
3312 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3313              if (sigma_theta(k,i).ne.0)
3314      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3315
3316 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3317 c                             rescore(k,i-2) !  right expression ?
3318 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3319           enddo
3320         endif
3321
3322         if (waga_d.gt.0.0d0) then
3323 c       open (ientin,file=tpl_k_sigma_d,status='old')
3324 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3325 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3326 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3327 c    &                          sigma_d(k,i+nnt-1)
3328 c         enddo
3329 c1404   continue
3330
3331           do i = nnt,nct ! right? without parallel.
3332 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
3333 c         do i=loc_start,loc_end ! with FG parallel.
3334                if (itype(i).eq.10) cycle
3335                if (idomain(k,i).eq.0 ) then
3336                   sigma_d(k,i)=0.0
3337                   cycle
3338                endif
3339                xxtpl(k,i)=xxref(i)
3340                yytpl(k,i)=yyref(i)
3341                zztpl(k,i)=zzref(i)
3342 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3343 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3344 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3345 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
3346                sigma_d(k,i)=rescore3(k,i) !  right expression ?
3347                if (sigma_d(k,i).ne.0)
3348      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3349
3350 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
3351 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3352 c              read (ientin,*) sigma_d(k,i) ! 1st variant
3353           enddo
3354         endif
3355       enddo
3356 c
3357 c remove distance restraints not used in any model from the list
3358 c shift data in all arrays
3359 c
3360       if (waga_dist.ne.0.0d0) then
3361         ii=0
3362         liiflag=.true.
3363         do i=nnt,nct-2
3364          do j=i+2,nct
3365           ii=ii+1
3366           if (ii_in_use(ii).eq.0.and.liiflag) then
3367             liiflag=.false.
3368             iistart=ii
3369           endif
3370           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3371      &                   .not.liiflag.and.ii.eq.lim_odl) then
3372              if (ii.eq.lim_odl) then
3373               iishift=ii-iistart+1
3374              else
3375               iishift=ii-iistart
3376              endif
3377              liiflag=.true.
3378              do ki=iistart,lim_odl-iishift
3379               ires_homo(ki)=ires_homo(ki+iishift)
3380               jres_homo(ki)=jres_homo(ki+iishift)
3381               ii_in_use(ki)=ii_in_use(ki+iishift)
3382               do k=1,constr_homology
3383                odl(k,ki)=odl(k,ki+iishift)
3384                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3385                l_homo(k,ki)=l_homo(k,ki+iishift)
3386               enddo
3387              enddo
3388              ii=ii-iishift
3389              lim_odl=lim_odl-iishift
3390           endif
3391          enddo
3392         enddo
3393       endif
3394
3395       endif ! .not. klapaucjusz
3396
3397       if (constr_homology.gt.0) call homology_partition
3398       if (constr_homology.gt.0) call init_int_table
3399 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3400 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3401 c
3402 c Print restraints
3403 c
3404       if (.not.out_template_restr) return
3405 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3406       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3407        write (iout,*) "Distance restraints from templates"
3408        do ii=1,lim_odl
3409        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3410      &  ii,ires_homo(ii),jres_homo(ii),
3411      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3412      &  ki=1,constr_homology)
3413        enddo
3414        write (iout,*) "Dihedral angle restraints from templates"
3415        do i=nnt+3,nct
3416         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3417      &      (rad2deg*dih(ki,i),
3418      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3419        enddo
3420        write (iout,*) "Virtual-bond angle restraints from templates"
3421        do i=nnt+2,nct
3422         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3423      &      (rad2deg*thetatpl(ki,i),
3424      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3425        enddo
3426        write (iout,*) "SC restraints from templates"
3427        do i=nnt,nct
3428         write(iout,'(i5,100(4f8.2,4x))') i,
3429      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3430      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3431        enddo
3432       endif
3433 c -----------------------------------------------------------------
3434       return
3435       end
3436 c----------------------------------------------------------------------
3437 #ifdef WINIFL
3438       subroutine flush(iu)
3439       return
3440       end
3441 #endif
3442 #ifdef AIX
3443       subroutine flush(iu)
3444       call flush_(iu)
3445       return
3446       end
3447 #endif
3448 c------------------------------------------------------------------------------
3449       subroutine copy_to_tmp(source)
3450       include "DIMENSIONS"
3451       include "COMMON.IOUNITS"
3452       character*(*) source
3453       character* 256 tmpfile
3454       integer ilen
3455       external ilen
3456       logical ex
3457       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3458       inquire(file=tmpfile,exist=ex)
3459       if (ex) then
3460         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3461      &   " to temporary directory..."
3462         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3463         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3464       endif
3465       return
3466       end
3467 c------------------------------------------------------------------------------
3468       subroutine move_from_tmp(source)
3469       include "DIMENSIONS"
3470       include "COMMON.IOUNITS"
3471       character*(*) source
3472       integer ilen
3473       external ilen
3474       write (*,*) "Moving ",source(:ilen(source)),
3475      & " from temporary directory to working directory"
3476       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3477       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3478       return
3479       end
3480 c------------------------------------------------------------------------------
3481       subroutine random_init(seed)
3482 C
3483 C Initialize random number generator
3484 C
3485       implicit real*8 (a-h,o-z)
3486       include 'DIMENSIONS'
3487 #ifdef MPI
3488       include 'mpif.h'
3489       logical OKRandom, prng_restart
3490       real*8  r1
3491       integer iseed_array(4)
3492 #endif
3493       include 'COMMON.IOUNITS'
3494       include 'COMMON.TIME1'
3495       include 'COMMON.CONTROL'
3496       include 'COMMON.SETUP'
3497       iseed=-dint(dabs(seed))
3498       if (iseed.eq.0) then
3499         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3500      &    'Random seed undefined. The program will stop.'
3501         write (*,'(/80(1h*)/20x,a/80(1h*))')
3502      &    'Random seed undefined. The program will stop.'
3503 #ifdef MPI
3504         call mpi_finalize(mpi_comm_world,ierr)
3505 #endif
3506         stop 'Bad random seed.'
3507       endif
3508 #ifdef MPI
3509       if (fg_rank.eq.0) then
3510       seed=seed*(me+1)+1
3511 #ifdef AMD64
3512       if(me.eq.king)
3513      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3514       OKRandom = prng_restart(me,iseed)
3515 #else
3516       do i=1,4
3517        tmp=65536.0d0**(4-i)
3518        iseed_array(i) = dint(seed/tmp)
3519        seed=seed-iseed_array(i)*tmp
3520       enddo
3521       if(me.eq.king)
3522      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3523      &                 (iseed_array(i),i=1,4)
3524       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3525      &                 (iseed_array(i),i=1,4)
3526       OKRandom = prng_restart(me,iseed_array)
3527 #endif
3528       if (OKRandom) then
3529 c        r1 = prng_next(me)
3530         r1=ran_number(0.0D0,1.0D0)
3531         if(me.eq.king)
3532      &   write (iout,*) 'ran_num',r1
3533         if (r1.lt.0.0d0) OKRandom=.false.
3534       endif
3535       if (.not.OKRandom) then
3536         write (iout,*) 'PRNG IS NOT WORKING!!!'
3537         print *,'PRNG IS NOT WORKING!!!'
3538         if (me.eq.0) then
3539          call flush(iout)
3540          call mpi_abort(mpi_comm_world,error_msg,ierr)
3541          stop
3542         else
3543          write (iout,*) 'too many processors for parallel prng'
3544          write (*,*) 'too many processors for parallel prng'
3545          call flush(iout)
3546          stop
3547         endif
3548       endif
3549       endif
3550 #else
3551       call vrndst(iseed)
3552       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3553 #endif
3554       return
3555       end
3556 c----------------------------------------------------------------------
3557       subroutine read_klapaucjusz
3558
3559       include 'DIMENSIONS'
3560 #ifdef MPI
3561       include 'mpif.h'
3562 #endif
3563       include 'COMMON.SETUP'
3564       include 'COMMON.CONTROL'
3565       include 'COMMON.CHAIN'
3566       include 'COMMON.IOUNITS'
3567       include 'COMMON.HOMOLOGY'
3568       include 'COMMON.GEO'
3569       include 'COMMON.INTERACT'
3570       include 'COMMON.NAMES'
3571       character*256 fragfile
3572       integer ninclust(maxclust),inclust(max_template,maxclust),
3573      &  nresclust(maxclust),iresclust(maxres,maxclust)
3574
3575       character*2 kic2
3576       character*24 model_ki_dist, model_ki_angle
3577       character*500 controlcard
3578       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3579       logical lprn /.true./
3580       integer ilen
3581       external ilen
3582       logical liiflag
3583 c
3584 c
3585       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3586       double precision, dimension (max_template,maxres) :: rescore
3587       double precision, dimension (max_template,maxres) :: rescore2
3588       character*24 pdbfile,tpl_k_rescore
3589
3590 c
3591 c For new homol impl
3592 c
3593       include 'COMMON.VAR'
3594 c
3595       call getenv("FRAGFILE",fragfile)
3596       open(ientin,file=fragfile,status="old",err=10)
3597       read(ientin,*) constr_homology,nclust
3598       l_homo = .false.
3599       sigma_theta=0.0
3600       sigma_d=0.0
3601       sigma_dih=0.0
3602 c Read pdb files
3603       do k=1,constr_homology
3604         read(ientin,'(a)') pdbfile
3605         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3606      &  pdbfile(:ilen(pdbfile))
3607         open(ipdbin,file=pdbfile,status='old',err=33)
3608         goto 34
3609   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3610      &  pdbfile(:ilen(pdbfile))
3611         stop
3612   34    continue
3613         unres_pdb=.false.
3614         call readpdb_template(k)
3615         do i=1,nres
3616           rescore(k,i)=0.2d0
3617           rescore2(k,i)=1.0d0
3618         enddo
3619       enddo
3620 c Read clusters
3621       do i=1,nclust
3622         read(ientin,*) ninclust(i),nresclust(i)
3623         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3624         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3625       enddo
3626 c
3627 c Loop over clusters
3628 c
3629       do l=1,nclust
3630         do ll = 1,ninclust(l)
3631
3632         k = inclust(ll,l)
3633         do i=1,nres
3634           idomain(k,i)=0
3635         enddo
3636         do i=1,nresclust(l)
3637           if (nnt.gt.1)  then
3638             idomain(k,iresclust(i,l)+1) = 1
3639           else
3640             idomain(k,iresclust(i,l)) = 1
3641           endif
3642         enddo
3643 c
3644 c     Distance restraints
3645 c
3646 c          ... --> odl(k,ii)
3647 C Copy the coordinates from reference coordinates (?)
3648         do i=1,2*nres
3649           do j=1,3
3650             c(j,i)=chomo(j,i,k)
3651 c           write (iout,*) "c(",j,i,") =",c(j,i)
3652           enddo
3653         enddo
3654         call int_from_cart(.true.,.false.)
3655         call sc_loc_geom(.false.)
3656         do i=1,nres
3657           thetaref(i)=theta(i)
3658           phiref(i)=phi(i)
3659         enddo
3660         if (waga_dist.ne.0.0d0) then
3661           ii=0
3662           do i = nnt,nct-2
3663             do j=i+2,nct
3664
3665               x12=c(1,i)-c(1,j)
3666               y12=c(2,i)-c(2,j)
3667               z12=c(3,i)-c(3,j)
3668               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3669 c              write (iout,*) k,i,j,distal,dist2_cut
3670
3671             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3672      &            .and. distal.le.dist2_cut ) then
3673
3674               ii=ii+1
3675               ii_in_use(ii)=1
3676               l_homo(k,ii)=.true.
3677
3678 c             write (iout,*) "k",k
3679 c             write (iout,*) "i",i," j",j," constr_homology",
3680 c    &                       constr_homology
3681               ires_homo(ii)=i
3682               jres_homo(ii)=j
3683               odl(k,ii)=distal
3684               if (read2sigma) then
3685                 sigma_odl(k,ii)=0
3686                 do ik=i,j
3687                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3688                 enddo
3689                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3690                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3691      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3692               else
3693                 if (odl(k,ii).le.dist_cut) then
3694                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3695                 else
3696 #ifdef OLDSIGMA
3697                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3698      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3699 #else
3700                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3701      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3702 #endif
3703                 endif
3704               endif
3705               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3706             else
3707               ii=ii+1
3708 c              l_homo(k,ii)=.false.
3709             endif
3710             enddo
3711           enddo
3712         lim_odl=ii
3713         endif
3714 c
3715 c     Theta, dihedral and SC retraints
3716 c
3717         if (waga_angle.gt.0.0d0) then
3718           do i = nnt+3,nct
3719             if (idomain(k,i).eq.0) then
3720 c               sigma_dih(k,i)=0.0
3721                cycle
3722             endif
3723             dih(k,i)=phiref(i)
3724             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3725      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3726 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3727 c     &       " sigma_dihed",sigma_dih(k,i)
3728             if (sigma_dih(k,i).ne.0)
3729      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3730           enddo
3731           lim_dih=nct-nnt-2
3732         endif
3733
3734         if (waga_theta.gt.0.0d0) then
3735           do i = nnt+2,nct
3736              if (idomain(k,i).eq.0) then
3737 c              sigma_theta(k,i)=0.0
3738               cycle
3739              endif
3740              thetatpl(k,i)=thetaref(i)
3741              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3742      &                        rescore(k,i-2))/3.0
3743              if (sigma_theta(k,i).ne.0)
3744      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3745           enddo
3746         endif
3747
3748         if (waga_d.gt.0.0d0) then
3749           do i = nnt,nct
3750                if (itype(i).eq.10) cycle
3751                if (idomain(k,i).eq.0 ) then
3752 c                  sigma_d(k,i)=0.0
3753                   cycle
3754                endif
3755                xxtpl(k,i)=xxref(i)
3756                yytpl(k,i)=yyref(i)
3757                zztpl(k,i)=zzref(i)
3758                sigma_d(k,i)=rescore(k,i)
3759                if (sigma_d(k,i).ne.0)
3760      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3761                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3762           enddo
3763         endif
3764       enddo ! l
3765       enddo ! ll
3766 c
3767 c remove distance restraints not used in any model from the list
3768 c shift data in all arrays
3769 c
3770       if (waga_dist.ne.0.0d0) then
3771         ii=0
3772         liiflag=.true.
3773         do i=nnt,nct-2
3774          do j=i+2,nct
3775           ii=ii+1
3776           if (ii_in_use(ii).eq.0.and.liiflag) then
3777             liiflag=.false.
3778             iistart=ii
3779           endif
3780           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3781      &                   .not.liiflag.and.ii.eq.lim_odl) then
3782              if (ii.eq.lim_odl) then
3783               iishift=ii-iistart+1
3784              else
3785               iishift=ii-iistart
3786              endif
3787              liiflag=.true.
3788              do ki=iistart,lim_odl-iishift
3789               ires_homo(ki)=ires_homo(ki+iishift)
3790               jres_homo(ki)=jres_homo(ki+iishift)
3791               ii_in_use(ki)=ii_in_use(ki+iishift)
3792               do k=1,constr_homology
3793                odl(k,ki)=odl(k,ki+iishift)
3794                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3795                l_homo(k,ki)=l_homo(k,ki+iishift)
3796               enddo
3797              enddo
3798              ii=ii-iishift
3799              lim_odl=lim_odl-iishift
3800           endif
3801          enddo
3802         enddo
3803       endif
3804
3805       return
3806    10 stop "Error infragment file"
3807       end