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