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