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