source/unres/src-HCD-5D/readrtns_CSA.F

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