2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.CHAIN'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 csa if (modecalc.eq.8) then
33 csa inquire (file="fort.40",exist=file_exist)
34 csa if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
46 & "The following",nhpb-nss,
47 & " distance restraints have been imposed:",
48 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
51 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
52 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
59 c print *,"Processor",myrank," leaves READRTNS"
62 C-------------------------------------------------------------------------------
63 subroutine read_control
67 implicit real*8 (a-h,o-z)
71 logical OKRandom, prng_restart
74 include 'COMMON.IOUNITS'
75 include 'COMMON.TIME1'
76 include 'COMMON.THREAD'
77 include 'COMMON.SBRIDGE'
78 include 'COMMON.CONTROL'
81 include 'COMMON.HEADER'
82 csa include 'COMMON.CSA'
83 include 'COMMON.CHAIN'
86 include 'COMMON.FFIELD'
87 include 'COMMON.SETUP'
88 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
89 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
91 character*320 controlcard
96 read (INP,'(a)') titel
97 call card_concat(controlcard)
98 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
99 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
100 call reada(controlcard,'SEED',seed,0.0D0)
101 call random_init(seed)
102 C Set up the time limit (caution! The time must be input in minutes!)
103 read_cart=index(controlcard,'READ_CART').gt.0
104 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
105 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
106 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
107 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
108 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
109 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
110 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
111 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
112 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
113 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
114 call reada(controlcard,'DRMS',drms,0.1D0)
115 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
116 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
117 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
118 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
119 write (iout,'(a,f10.1)')'DRMS = ',drms
120 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
121 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
123 call readi(controlcard,'NZ_START',nz_start,0)
124 call readi(controlcard,'NZ_END',nz_end,0)
125 call readi(controlcard,'IZ_SC',iz_sc,0)
127 safety = 60.0d0*safety
130 call reada(controlcard,"T_BATH",t_bath,300.0d0)
131 minim=(index(controlcard,'MINIMIZE').gt.0)
132 dccart=(index(controlcard,'CART').gt.0)
133 overlapsc=(index(controlcard,'OVERLAP').gt.0)
134 overlapsc=.not.overlapsc
135 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
136 searchsc=.not.searchsc
137 sideadd=(index(controlcard,'SIDEADD').gt.0)
138 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
139 outpdb=(index(controlcard,'PDBOUT').gt.0)
140 outx=(index(controlcard,'XOUT').gt.0)
141 outmol2=(index(controlcard,'MOL2OUT').gt.0)
142 pdbref=(index(controlcard,'PDBREF').gt.0)
143 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
144 indpdb=index(controlcard,'PDBSTART')
145 extconf=(index(controlcard,'EXTCONF').gt.0)
146 call readi(controlcard,'IPRINT',iprint,0)
147 call readi(controlcard,'MAXGEN',maxgen,10000)
148 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
149 call readi(controlcard,"KDIAG",kdiag,0)
150 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
151 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
152 & write (iout,*) "RESCALE_MODE",rescale_mode
153 split_ene=index(controlcard,'SPLIT_ENE').gt.0
154 if (index(controlcard,'REGULAR').gt.0.0D0) then
155 call reada(controlcard,'WEIDIS',weidis,0.1D0)
159 call reada(controlcard,"CHECKGRAD_INC",checkgrad_inc,1.0d-4)
160 if (index(controlcard,'CHECKGRAD').gt.0) then
162 if (index(controlcard,'CART').gt.0) then
164 elseif (index(controlcard,'CARINT').gt.0) then
169 elseif (index(controlcard,'THREAD').gt.0) then
171 call readi(controlcard,'THREAD',nthread,0)
172 if (nthread.gt.0) then
173 call reada(controlcard,'WEIDIS',weidis,0.1D0)
176 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
177 stop 'Error termination in Read_Control.'
179 else if (index(controlcard,'MCMA').gt.0) then
181 else if (index(controlcard,'MCEE').gt.0) then
183 else if (index(controlcard,'MULTCONF').gt.0) then
185 else if (index(controlcard,'MAP').gt.0) then
187 call readi(controlcard,'MAP',nmap,0)
188 else if (index(controlcard,'CSA').gt.0) then
189 write(*,*) "CSA not supported in this version"
192 crc else if (index(controlcard,'ZSCORE').gt.0) then
194 crc ZSCORE is rm from UNRES, modecalc=9 is available
197 cfcm else if (index(controlcard,'MCMF').gt.0) then
199 else if (index(controlcard,'SOFTREG').gt.0) then
201 else if (index(controlcard,'CHECK_BOND').gt.0) then
203 else if (index(controlcard,'TEST').gt.0) then
205 else if (index(controlcard,'MD').gt.0) then
207 else if (index(controlcard,'RE ').gt.0) then
211 lmuca=index(controlcard,'MUCA').gt.0
212 call readi(controlcard,'MUCADYN',mucadyn,0)
213 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
214 if (lmuca .and. (me.eq.king .or. .not.out1file ))
216 write (iout,*) 'MUCADYN=',mucadyn
217 write (iout,*) 'MUCASMOOTH=',muca_smooth
220 iscode=index(controlcard,'ONE_LETTER')
221 indphi=index(controlcard,'PHI')
222 indback=index(controlcard,'BACK')
223 iranconf=index(controlcard,'RAND_CONF')
224 i2ndstr=index(controlcard,'USE_SEC_PRED')
225 gradout=index(controlcard,'GRADOUT').gt.0
226 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
228 if(me.eq.king.or..not.out1file)
229 & write (iout,'(2a)') diagmeth(kdiag),
230 & ' routine used to diagonalize matrices.'
233 c--------------------------------------------------------------------------
234 subroutine read_REMDpar
238 implicit real*8 (a-h,o-z)
240 include 'COMMON.IOUNITS'
241 include 'COMMON.TIME1'
244 include 'COMMON.LANGEVIN'
246 include 'COMMON.LANGEVIN.lang0'
248 include 'COMMON.INTERACT'
249 include 'COMMON.NAMES'
251 include 'COMMON.REMD'
252 include 'COMMON.CONTROL'
253 include 'COMMON.SETUP'
255 character*320 controlcard
256 character*3200 controlcard1
257 integer iremd_m_total
259 if(me.eq.king.or..not.out1file)
260 & write (iout,*) "REMD setup"
262 call card_concat(controlcard)
263 call readi(controlcard,"NREP",nrep,3)
264 call readi(controlcard,"NSTEX",nstex,1000)
265 call reada(controlcard,"RETMIN",retmin,10.0d0)
266 call reada(controlcard,"RETMAX",retmax,1000.0d0)
267 mremdsync=(index(controlcard,'SYNC').gt.0)
268 call readi(controlcard,"NSYN",i_sync_step,100)
269 restart1file=(index(controlcard,'REST1FILE').gt.0)
270 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
271 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
272 if(max_cache_traj_use.gt.max_cache_traj)
273 & max_cache_traj_use=max_cache_traj
274 if(me.eq.king.or..not.out1file) then
275 cd if (traj1file) then
276 crc caching is in testing - NTWX is not ignored
277 cd write (iout,*) "NTWX value is ignored"
278 cd write (iout,*) " trajectory is stored to one file by master"
279 cd write (iout,*) " before exchange at NSTEX intervals"
281 write (iout,*) "NREP= ",nrep
282 write (iout,*) "NSTEX= ",nstex
283 write (iout,*) "SYNC= ",mremdsync
284 write (iout,*) "NSYN= ",i_sync_step
285 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
288 t_exchange_only=(index(controlcard,'TONLY').gt.0)
289 call readi(controlcard,"HREMD",hremd,0)
290 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
291 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
293 if(usampl.and.hremd.gt.0) then
295 & "========== ERROR: USAMPL and HREMD cannot be used together"
297 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
304 if (index(controlcard,'TLIST').gt.0) then
306 call card_concat(controlcard1)
307 read(controlcard1,*) (remd_t(i),i=1,nrep)
308 if(me.eq.king.or..not.out1file)
309 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
312 if (index(controlcard,'MLIST').gt.0) then
314 call card_concat(controlcard1)
315 read(controlcard1,*) (remd_m(i),i=1,nrep)
316 if(me.eq.king.or..not.out1file) then
317 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
320 iremd_m_total=iremd_m_total+remd_m(i)
323 write (iout,*) 'Total number of replicas ',
324 & iremd_m_total*hremd
326 write (iout,*) 'Total number of replicas ',iremd_m_total
330 if(me.eq.king.or..not.out1file)
331 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
334 c--------------------------------------------------------------------------
335 subroutine read_MDpar
339 implicit real*8 (a-h,o-z)
341 include 'COMMON.IOUNITS'
342 include 'COMMON.TIME1'
345 include 'COMMON.LANGEVIN'
347 include 'COMMON.LANGEVIN.lang0'
349 include 'COMMON.INTERACT'
350 include 'COMMON.NAMES'
352 include 'COMMON.SETUP'
353 include 'COMMON.CONTROL'
354 include 'COMMON.SPLITELE'
356 character*320 controlcard
358 call card_concat(controlcard)
359 call readi(controlcard,"NSTEP",n_timestep,1000000)
360 call readi(controlcard,"NTWE",ntwe,100)
361 call readi(controlcard,"NTWX",ntwx,1000)
362 call reada(controlcard,"DT",d_time,1.0d-1)
363 call reada(controlcard,"DVMAX",dvmax,2.0d1)
364 call reada(controlcard,"DAMAX",damax,1.0d1)
365 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
366 call readi(controlcard,"LANG",lang,0)
367 RESPA = index(controlcard,"RESPA") .gt. 0
368 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
369 ntime_split0=ntime_split
370 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
371 ntime_split0=ntime_split
372 call reada(controlcard,"R_CUT",r_cut,2.0d0)
373 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
374 rest = index(controlcard,"REST").gt.0
375 tbf = index(controlcard,"TBF").gt.0
376 call readi(controlcard,"HMC",hmc,0)
377 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
378 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
379 tnh = index(controlcard,"NOSEHOOVER96").gt.0
380 if (RESPA.and.tnh)then
381 xiresp = index(controlcard,"XIRESP").gt.0
383 call reada(controlcard,"Q_NP",Q_np,0.1d0)
384 usampl = index(controlcard,"USAMPL").gt.0
385 mdpdb = index(controlcard,"MDPDB").gt.0
386 call reada(controlcard,"T_BATH",t_bath,300.0d0)
387 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
388 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
389 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
390 if (count_reset_moment.eq.0) count_reset_moment=1000000000
391 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
392 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
393 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
394 if (count_reset_vel.eq.0) count_reset_vel=1000000000
395 large = index(controlcard,"LARGE").gt.0
396 print_compon = index(controlcard,"PRINT_COMPON").gt.0
397 rattle = index(controlcard,"RATTLE").gt.0
398 preminim = index(controlcard,"PREMINIM").gt.0
400 dccart=(index(controlcard,'CART').gt.0)
403 c if performing umbrella sampling, fragments constrained are read from the fragment file
409 if(me.eq.king.or..not.out1file) then
411 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
413 write (iout,'(a)') "The units are:"
414 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
415 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
416 & " acceleration: angstrom/(48.9 fs)**2"
417 write (iout,'(a)') "energy: kcal/mol, temperature: K"
419 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
420 write (iout,'(a60,f10.5,a)')
421 & "Initial time step of numerical integration:",d_time,
423 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
425 write (iout,'(2a,i4,a)')
426 & "A-MTS algorithm used; initial time step for fast-varying",
427 & " short-range forces split into",ntime_split," steps."
428 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
429 & r_cut," lambda",rlamb
431 write (iout,'(2a,f10.5)')
432 & "Maximum acceleration threshold to reduce the time step",
433 & "/increase split number:",damax
434 write (iout,'(2a,f10.5)')
435 & "Maximum predicted energy drift to reduce the timestep",
436 & "/increase split number:",edriftmax
437 write (iout,'(a60,f10.5)')
438 & "Maximum velocity threshold to reduce velocities:",dvmax
439 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
440 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
441 if (rattle) write (iout,'(a60)')
442 & "Rattle algorithm used to constrain the virtual bonds"
443 if (preminim .or. iranconf.gt.0) then
445 & "Initial structure will be energy-minimized"
450 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
451 call reada(controlcard,"RWAT",rwat,1.4d0)
452 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
453 surfarea=index(controlcard,"SURFAREA").gt.0
454 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
455 if(me.eq.king.or..not.out1file)then
456 write (iout,'(/a,$)') "Langevin dynamics calculation"
459 & " with direct integration of Langevin equations"
460 else if (lang.eq.2) then
461 write (iout,'(a/)') " with TINKER stochasic MD integrator"
462 else if (lang.eq.3) then
463 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
464 else if (lang.eq.4) then
465 write (iout,'(a/)') " in overdamped mode"
467 write (iout,'(//a,i5)')
468 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
471 write (iout,'(a60,f10.5)') "Temperature:",t_bath
472 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
473 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
474 write (iout,'(a60,f10.5)')
475 & "Scaling factor of the friction forces:",scal_fric
476 if (surfarea) write (iout,'(2a,i10,a)')
477 & "Friction coefficients will be scaled by solvent-accessible",
478 & " surface area every",reset_fricmat," steps."
480 c Calculate friction coefficients and bounds of stochastic forces
481 eta=6*pi*cPoise*etawat
482 if(me.eq.king.or..not.out1file)
483 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
485 gamp=scal_fric*(pstok+rwat)*eta
486 stdfp=dsqrt(2*Rb*t_bath/d_time)
488 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
489 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
491 if(me.eq.king.or..not.out1file)then
492 write (iout,'(/2a/)')
493 & "Radii of site types and friction coefficients and std's of",
494 & " stochastic forces of fully exposed sites"
495 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
497 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
498 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
502 if(me.eq.king.or..not.out1file)then
503 write (iout,'(a)') "Berendsen bath calculation"
504 write (iout,'(a60,f10.5)') "Temperature:",t_bath
505 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
507 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
508 & count_reset_moment," steps"
510 & write (iout,'(a,i10,a)')
511 & "Velocities will be reset at random every",count_reset_vel,
514 else if (tnp .or. tnp1 .or. tnh) then
515 if (tnp .or. tnp1) then
516 write (iout,'(a)') "Nose-Poincare bath calculation"
517 if (tnp) write (iout,'(a)')
518 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
519 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
521 write (iout,'(a)') "Nose-Hoover bath calculation"
522 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
532 WDTI(i) = 1.0*d_time/nresn
539 write (iout,'(a)') "NVT-XI-RESPA algorithm"
541 write (iout,'(a)') "NVT-XO-RESPA algorithm"
544 WDTIi(i) = 1.0*d_time/nresn/ntime_split
552 write (iout,'(a60,f10.5)') "Temperature:",t_bath
553 write (iout,'(a60,f10.5)') "Q =",Q_np
555 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
556 & count_reset_moment," steps"
558 & write (iout,'(a,i10,a)')
559 & "Velocities will be reset at random every",count_reset_vel,
562 else if (hmc.gt.0) then
563 write (iout,'(a)') "Hybrid Monte Carlo calculation"
564 write (iout,'(a60,f10.5)') "Temperature:",t_bath
565 write (iout,'(a60,i10)')
566 & "Number of MD steps between Metropolis tests:",hmc
569 if(me.eq.king.or..not.out1file)
570 & write (iout,'(a31)') "Microcanonical mode calculation"
572 if(me.eq.king.or..not.out1file)then
573 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
575 write(iout,*) "MD running with constraints."
576 write(iout,*) "Equilibration time ", eq_time, " mtus."
577 write(iout,*) "Constraining ", nfrag," fragments."
578 write(iout,*) "Length of each fragment, weight and q0:"
580 write (iout,*) "Set of restraints #",iset
582 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
583 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
585 write(iout,*) "constraints between ", npair, "fragments."
586 write(iout,*) "constraint pairs, weights and q0:"
588 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
589 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
591 write(iout,*) "angle constraints within ", nfrag_back,
592 & "backbone fragments."
593 write(iout,*) "fragment, weights:"
595 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
596 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
597 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
600 iset=mod(kolor,nset)+1
603 if(me.eq.king.or..not.out1file)
604 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
607 c------------------------------------------------------------------------------
610 C Read molecular data.
612 implicit real*8 (a-h,o-z)
618 include 'COMMON.IOUNITS'
621 include 'COMMON.INTERACT'
622 include 'COMMON.LOCAL'
623 include 'COMMON.NAMES'
624 include 'COMMON.CHAIN'
625 include 'COMMON.FFIELD'
626 include 'COMMON.SBRIDGE'
627 include 'COMMON.HEADER'
628 include 'COMMON.CONTROL'
629 include 'COMMON.DBASE'
630 include 'COMMON.THREAD'
631 include 'COMMON.CONTACTS'
632 include 'COMMON.TORCNSTR'
633 include 'COMMON.TIME1'
634 include 'COMMON.BOUNDS'
636 include 'COMMON.REMD'
637 include 'COMMON.SETUP'
638 character*4 sequence(maxres)
640 double precision x(maxvar)
641 character*256 pdbfile
642 character*320 weightcard
643 character*80 weightcard_t,ucase
644 dimension itype_pdb(maxres)
645 common /pizda/ itype_pdb
646 logical seq_comp,fail
647 double precision energia(0:n_ene)
654 C Read weights of the subsequent energy terms.
667 if(me.eq.king.or..not.out1file) then
668 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
669 write (iout,*) 'Current weights for processor ',
670 & me,' set ',i2set(me)
674 call card_concat(weightcard)
675 call reada(weightcard,'WLONG',wlong,1.0D0)
676 call reada(weightcard,'WSC',wsc,wlong)
677 call reada(weightcard,'WSCP',wscp,wlong)
678 call reada(weightcard,'WELEC',welec,1.0D0)
679 call reada(weightcard,'WVDWPP',wvdwpp,welec)
680 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
681 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
682 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
683 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
684 call reada(weightcard,'WTURN3',wturn3,1.0D0)
685 call reada(weightcard,'WTURN4',wturn4,1.0D0)
686 call reada(weightcard,'WTURN6',wturn6,1.0D0)
687 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
688 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
689 call reada(weightcard,'WBOND',wbond,1.0D0)
690 call reada(weightcard,'WTOR',wtor,1.0D0)
691 call reada(weightcard,'WTORD',wtor_d,1.0D0)
692 call reada(weightcard,'WANG',wang,1.0D0)
693 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
694 call reada(weightcard,'SCAL14',scal14,0.4D0)
695 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
696 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
697 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
698 call reada(weightcard,'TEMP0',temp0,300.0d0)
699 if (index(weightcard,'SOFT').gt.0) ipot=6
700 C 12/1/95 Added weight for the multi-body term WCORR
701 call reada(weightcard,'WCORRH',wcorr,1.0D0)
702 if (wcorr4.gt.0.0d0) wcorr=wcorr4
710 hweights(i,7)=wel_loc
713 hweights(i,10)=wturn6
715 hweights(i,12)=wscloc
717 hweights(i,14)=wtor_d
718 hweights(i,15)=wstrain
719 hweights(i,16)=wvdwpp
721 hweights(i,18)=scal14
722 hweights(i,21)=wsccor
727 weights(i)=hweights(i2set(me),i)
751 call card_concat(weightcard)
752 call reada(weightcard,'WLONG',wlong,1.0D0)
753 call reada(weightcard,'WSC',wsc,wlong)
754 call reada(weightcard,'WSCP',wscp,wlong)
755 call reada(weightcard,'WELEC',welec,1.0D0)
756 call reada(weightcard,'WVDWPP',wvdwpp,welec)
757 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
758 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
759 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
760 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
761 call reada(weightcard,'WTURN3',wturn3,1.0D0)
762 call reada(weightcard,'WTURN4',wturn4,1.0D0)
763 call reada(weightcard,'WTURN6',wturn6,1.0D0)
764 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
765 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
766 call reada(weightcard,'WBOND',wbond,1.0D0)
767 call reada(weightcard,'WTOR',wtor,1.0D0)
768 call reada(weightcard,'WTORD',wtor_d,1.0D0)
769 call reada(weightcard,'WANG',wang,1.0D0)
770 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
771 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
772 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
773 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
774 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
775 call reada(weightcard,'SCAL14',scal14,0.4D0)
776 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
777 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
778 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
779 call reada(weightcard,'TEMP0',temp0,300.0d0)
780 if (index(weightcard,'SOFT').gt.0) ipot=6
781 C 12/1/95 Added weight for the multi-body term WCORR
782 call reada(weightcard,'WCORRH',wcorr,1.0D0)
783 if (wcorr4.gt.0.0d0) wcorr=wcorr4
804 weights(25)=wdfa_dist
807 weights(28)=wdfa_beta
809 if(me.eq.king.or..not.out1file)
810 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
811 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
813 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
814 10 format (/'Energy-term weights (unscaled):'//
815 & 'WSCC= ',f10.6,' (SC-SC)'/
816 & 'WSCP= ',f10.6,' (SC-p)'/
817 & 'WELEC= ',f10.6,' (p-p electr)'/
818 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
819 & 'WBOND= ',f10.6,' (stretching)'/
820 & 'WANG= ',f10.6,' (bending)'/
821 & 'WSCLOC= ',f10.6,' (SC local)'/
822 & 'WTOR= ',f10.6,' (torsional)'/
823 & 'WTORD= ',f10.6,' (double torsional)'/
824 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
825 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
826 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
827 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
828 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
829 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
830 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
831 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
832 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
833 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
834 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
835 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
836 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
837 if(me.eq.king.or..not.out1file)then
838 if (wcorr4.gt.0.0d0) then
839 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
840 & 'between contact pairs of peptide groups'
841 write (iout,'(2(a,f5.3/))')
842 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
843 & 'Range of quenching the correlation terms:',2*delt_corr
844 else if (wcorr.gt.0.0d0) then
845 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
846 & 'between contact pairs of peptide groups'
848 write (iout,'(a,f8.3)')
849 & 'Scaling factor of 1,4 SC-p interactions:',scal14
850 write (iout,'(a,f8.3)')
851 & 'General scaling factor of SC-p interactions:',scalscp
853 r0_corr=cutoff_corr-delt_corr
855 aad(i,1)=scalscp*aad(i,1)
856 aad(i,2)=scalscp*aad(i,2)
857 bad(i,1)=scalscp*bad(i,1)
858 bad(i,2)=scalscp*bad(i,2)
860 call rescale_weights(t_bath)
861 if(me.eq.king.or..not.out1file)
862 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
863 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
865 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
866 22 format (/'Energy-term weights (scaled):'//
867 & 'WSCC= ',f10.6,' (SC-SC)'/
868 & 'WSCP= ',f10.6,' (SC-p)'/
869 & 'WELEC= ',f10.6,' (p-p electr)'/
870 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
871 & 'WBOND= ',f10.6,' (stretching)'/
872 & 'WANG= ',f10.6,' (bending)'/
873 & 'WSCLOC= ',f10.6,' (SC local)'/
874 & 'WTOR= ',f10.6,' (torsional)'/
875 & 'WTORD= ',f10.6,' (double torsional)'/
876 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
877 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
878 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
879 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
880 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
881 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
882 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
883 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
884 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
885 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
886 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
887 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
888 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
889 if(me.eq.king.or..not.out1file)
890 & write (iout,*) "Reference temperature for weights calculation:",
892 call reada(weightcard,"D0CM",d0cm,3.78d0)
893 call reada(weightcard,"AKCM",akcm,15.1d0)
894 call reada(weightcard,"AKTH",akth,11.0d0)
895 call reada(weightcard,"AKCT",akct,12.0d0)
896 call reada(weightcard,"V1SS",v1ss,-1.08d0)
897 call reada(weightcard,"V2SS",v2ss,7.61d0)
898 call reada(weightcard,"V3SS",v3ss,13.7d0)
899 call reada(weightcard,"EBR",ebr,-5.50D0)
900 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
902 dyn_ss_mask(i)=.false.
906 dyn_ssbond_ij(i,j)=1.0d300
909 call reada(weightcard,"HT",Ht,0.0D0)
911 ss_depth=ebr/wsc-0.25*eps(1,1)
912 Ht=Ht/wsc-0.25*eps(1,1)
913 akcm=akcm*wstrain/wsc
914 akth=akth*wstrain/wsc
915 akct=akct*wstrain/wsc
916 v1ss=v1ss*wstrain/wsc
917 v2ss=v2ss*wstrain/wsc
918 v3ss=v3ss*wstrain/wsc
920 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
923 if(me.eq.king.or..not.out1file) then
924 write (iout,*) "Parameters of the SS-bond potential:"
925 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
927 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
928 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
929 write (iout,*)" HT",Ht
930 print *,'indpdb=',indpdb,' pdbref=',pdbref
932 if (indpdb.gt.0 .or. pdbref) then
933 read(inp,'(a)') pdbfile
934 if(me.eq.king.or..not.out1file)
935 & write (iout,'(2a)') 'PDB data will be read from file ',
936 & pdbfile(:ilen(pdbfile))
937 open(ipdbin,file=pdbfile,status='old',err=33)
939 33 write (iout,'(a)') 'Error opening PDB file.'
942 c print *,'Begin reading pdb data'
951 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
952 & (crefjlee(j,i+nres),j=1,3)
955 c print *,'Finished reading pdb data'
956 if(me.eq.king.or..not.out1file)
957 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
958 & ' nstart_sup=',nstart_sup
960 itype_pdb(i)=itype(i)
964 nct=nstart_sup+nsup-1
965 call contact(.false.,ncont_ref,icont_ref,co)
968 C Following 2 lines for diagnostics; comment out if not needed
969 write (iout,*) "Before sideadd"
971 if(me.eq.king.or..not.out1file)
972 & write(iout,*)'Adding sidechains'
979 do while (fail.and.nsi.le.maxsi)
980 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
983 if(fail) write(iout,*)'Adding sidechain failed for res ',
984 & i,' after ',nsi,' trials'
987 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
990 C Following 2 lines for diagnostics; comment out if not needed
991 c write (iout,*) "After sideadd"
994 if (indpdb.eq.0) then
995 C Read sequence if not taken from the pdb file.
997 c print *,'nres=',nres
998 if (iscode.gt.0) then
999 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1001 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1003 C Convert sequence to numeric code
1005 itype(i)=rescode(i,sequence(i),iscode)
1007 if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
1009 & "Glycine is the first full residue, initial dummy deleted"
1015 if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
1017 & "Glycine is the last full residue, terminal dummy deleted"
1021 C Assign initial virtual bond lengths
1027 vbld(i+nres)=dsc(itype(i))
1028 vbld_inv(i+nres)=dsc_inv(itype(i))
1029 c write (iout,*) "i",i," itype",itype(i),
1030 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
1034 c print '(20i4)',(itype(i),i=1,nres)
1037 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
1039 if (itype(i).eq.21) then
1043 else if (itype(i+1).ne.20) then
1045 else if (itype(i).ne.20) then
1052 if(me.eq.king.or..not.out1file)then
1053 write (iout,*) "ITEL"
1055 write (iout,*) i,itype(i),itel(i)
1057 print *,'Call Read_Bridge.'
1060 C 8/13/98 Set limits to generating the dihedral angles
1065 read (inp,*) ndih_constr
1066 if (ndih_constr.gt.0) then
1068 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1069 if(me.eq.king.or..not.out1file)then
1071 & 'There are',ndih_constr,' constraints on phi angles.'
1073 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1077 phi0(i)=deg2rad*phi0(i)
1078 drange(i)=deg2rad*drange(i)
1080 if(me.eq.king.or..not.out1file)
1081 & write (iout,*) 'FTORS',ftors
1084 phibound(1,ii) = phi0(i)-drange(i)
1085 phibound(2,ii) = phi0(i)+drange(i)
1090 if (me.eq.king) then
1092 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1094 write (iout,'(a3,i5,2f10.1)')
1095 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1101 cd print *,'NNT=',NNT,' NCT=',NCT
1102 if (itype(1).eq.21) nnt=2
1103 if (itype(nres).eq.21) nct=nct-1
1105 C Bartek:READ init_vars
1106 C Initialize variables!
1107 C Juyong:READ read_info
1108 C READ fragment information!!
1109 C both routines should be in dfa.F file!!
1112 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1113 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1115 print*, 'init_dfa_vars finished!'
1117 print*, 'read_dfa_info finished!'
1122 if(me.eq.king.or..not.out1file)
1123 & write (iout,'(a,i3)') 'nsup=',nsup
1125 if (nsup.le.(nct-nnt+1)) then
1126 do i=0,nct-nnt+1-nsup
1127 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1133 & 'Error - sequences to be superposed do not match.'
1136 do i=0,nsup-(nct-nnt+1)
1137 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1139 nstart_sup=nstart_sup+i
1145 & 'Error - sequences to be superposed do not match.'
1148 if (nsup.eq.0) nsup=nct-nnt
1149 if (nstart_sup.eq.0) nstart_sup=nnt
1150 if (nstart_seq.eq.0) nstart_seq=nnt
1151 if(me.eq.king.or..not.out1file)
1152 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1153 & ' nstart_seq=',nstart_seq
1155 c--- Zscore rms -------
1156 if (nz_start.eq.0) nz_start=nnt
1157 if (nz_end.eq.0 .and. nsup.gt.0) then
1159 else if (nz_end.eq.0) then
1162 if(me.eq.king.or..not.out1file)then
1163 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1164 write (iout,*) 'IZ_SC=',iz_sc
1166 c----------------------
1169 if (.not.pdbref) then
1170 call read_angles(inp,*38)
1172 38 write (iout,'(a)') 'Error reading reference structure.'
1174 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1175 stop 'Error reading reference structure'
1179 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1188 call contact(.true.,ncont_ref,icont_ref,co)
1190 if(me.eq.king.or..not.out1file)
1191 & write (iout,*) 'Contact order:',co
1193 if(me.eq.king.or..not.out1file)
1194 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1197 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1199 if(me.eq.king.or..not.out1file)
1200 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1201 & icont_ref(1,i),' ',
1202 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1206 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1207 if (constr_dist.gt.0) then
1208 call read_dist_constr
1212 if (constr_homology.gt.0) then
1213 call read_constr_homology
1214 if (indpdb.gt.0 .or. pdbref) then
1217 c(j,i)=crefjlee(j,i)
1218 cref(j,i)=crefjlee(j,i)
1223 write (iout,*) "Array C"
1225 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1226 & (c(j,i+nres),j=1,3)
1228 write (iout,*) "Array Cref"
1230 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1231 & (cref(j,i+nres),j=1,3)
1234 call int_from_cart1(.false.)
1235 call sc_loc_geom(.false.)
1237 thetaref(i)=theta(i)
1242 dc(j,i)=c(j,i+1)-c(j,i)
1243 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1248 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1249 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1257 if (nhpb.gt.0) call hpb_partition
1258 c write (iout,*) "After read_dist_constr nhpb",nhpb
1260 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1261 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1262 & modecalc.ne.10) then
1263 C If input structure hasn't been supplied from the PDB file read or generate
1265 if (iranconf.eq.0 .and. .not. extconf) then
1266 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1267 & write (iout,'(a)') 'Initial geometry will be read in.'
1269 c read (inp,*) time,potE,uconst,t_bath,
1270 c & nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn)
1271 read(inp,'(8f10.5)',end=36,err=36)
1272 & ((c(l,k),l=1,3),k=1,nres),
1273 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1274 call int_from_cart1(.false.)
1277 dc(j,i)=c(j,i+1)-c(j,i)
1278 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1282 if (itype(i).ne.10) then
1284 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1285 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1291 call read_angles(inp,*36)
1294 36 write (iout,'(a)') 'Error reading angle file.'
1296 call mpi_finalize( MPI_COMM_WORLD,IERR )
1298 stop 'Error reading angle file.'
1300 else if (extconf) then
1301 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1302 & write (iout,'(a)') 'Extended chain initial geometry.'
1304 theta(i)=90d0*deg2rad
1307 phi(i)=180d0*deg2rad
1310 alph(i)=110d0*deg2rad
1313 omeg(i)=-120d0*deg2rad
1316 if(me.eq.king.or..not.out1file)
1317 & write (iout,'(a)') 'Random-generated initial geometry.'
1321 if (me.eq.king .or. fg_rank.eq.0 .and. (
1322 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1326 call gen_rand_conf(itmp,*30)
1328 30 write (iout,*) 'Failed to generate random conformation',
1329 & ', itrial=',itrial
1330 write (*,*) 'Processor:',me,
1331 & ' Failed to generate random conformation',
1341 write (iout,'(a,i3,a)') 'Processor:',me,
1342 & ' error in generating random conformation.'
1343 write (*,'(a,i3,a)') 'Processor:',me,
1344 & ' error in generating random conformation.'
1347 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1354 elseif (modecalc.eq.4) then
1355 read (inp,'(a)') intinname
1356 open (intin,file=intinname,status='old',err=333)
1357 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1358 & write (iout,'(a)') 'intinname',intinname
1359 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1361 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1363 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1365 stop 'Error opening angle file.'
1369 C Generate distance constraints, if the PDB structure is to be regularized.
1370 if (nthread.gt.0) then
1371 call read_threadbase
1373 write (iout,*) "READRTNS: Calling setup_var"
1376 if (me.eq.king .or. .not. out1file)
1378 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1379 write (iout,'(/a,i3,a)')
1380 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1381 write (iout,'(20i4)') (iss(i),i=1,ns)
1383 write(iout,*)"Running with dynamic disulfide-bond formation"
1385 write (iout,'(/a/)') 'Pre-formed links are:'
1391 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1392 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1398 if (ns.gt.0.and.dyn_ss) then
1402 forcon(i-nss)=forcon(i)
1409 dyn_ss_mask(iss(i))=.true.
1412 if (i2ndstr.gt.0) call secstrp2dihc
1413 c call geom_to_var(nvar,x)
1414 c call etotal(energia(0))
1415 c call enerprint(energia(0))
1416 c call briefout(0,etot)
1418 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1419 cd write (iout,'(a)') 'Variable list:'
1420 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1422 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1423 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1424 & 'Processor',myrank,': end reading molecular data.'
1428 c--------------------------------------------------------------------------
1429 logical function seq_comp(itypea,itypeb,length)
1431 integer length,itypea(length),itypeb(length)
1434 if (itypea(i).ne.itypeb(i)) then
1442 c-----------------------------------------------------------------------------
1443 subroutine read_bridge
1444 C Read information about disulfide bridges.
1445 implicit real*8 (a-h,o-z)
1446 include 'DIMENSIONS'
1450 include 'COMMON.IOUNITS'
1451 include 'COMMON.GEO'
1452 include 'COMMON.VAR'
1453 include 'COMMON.INTERACT'
1454 include 'COMMON.LOCAL'
1455 include 'COMMON.NAMES'
1456 include 'COMMON.CHAIN'
1457 include 'COMMON.FFIELD'
1458 include 'COMMON.SBRIDGE'
1459 include 'COMMON.HEADER'
1460 include 'COMMON.CONTROL'
1461 include 'COMMON.DBASE'
1462 include 'COMMON.THREAD'
1463 include 'COMMON.TIME1'
1464 include 'COMMON.SETUP'
1465 C Read bridging residues.
1466 read (inp,*) ns,(iss(i),i=1,ns)
1468 if(me.eq.king.or..not.out1file)
1469 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1470 C Check whether the specified bridging residues are cystines.
1472 if (itype(iss(i)).ne.1) then
1473 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1474 & 'Do you REALLY think that the residue ',
1475 & restyp(itype(iss(i))),i,
1476 & ' can form a disulfide bridge?!!!'
1477 write (*,'(2a,i3,a)')
1478 & 'Do you REALLY think that the residue ',
1479 & restyp(itype(iss(i))),i,
1480 & ' can form a disulfide bridge?!!!'
1482 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1487 C Read preformed bridges.
1489 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1491 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1494 C Check if the residues involved in bridges are in the specified list of
1495 C bridging residues.
1498 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1499 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1500 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1501 & ' contains residues present in other pairs.'
1502 write (*,'(a,i3,a)') 'Disulfide pair',i,
1503 & ' contains residues present in other pairs.'
1505 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1511 if (ihpb(i).eq.iss(j)) goto 10
1513 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1516 if (jhpb(i).eq.iss(j)) goto 20
1518 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1524 ihpb(i)=ihpb(i)+nres
1525 jhpb(i)=jhpb(i)+nres
1531 c----------------------------------------------------------------------------
1532 subroutine read_x(kanal,*)
1533 implicit real*8 (a-h,o-z)
1534 include 'DIMENSIONS'
1535 include 'COMMON.GEO'
1536 include 'COMMON.VAR'
1537 include 'COMMON.CHAIN'
1538 include 'COMMON.IOUNITS'
1539 include 'COMMON.CONTROL'
1540 include 'COMMON.LOCAL'
1541 include 'COMMON.INTERACT'
1542 c Read coordinates from input
1544 read(kanal,'(8f10.5)',end=10,err=10)
1545 & ((c(l,k),l=1,3),k=1,nres),
1546 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1549 c(j,2*nres)=c(j,nres)
1551 call int_from_cart1(.false.)
1554 dc(j,i)=c(j,i+1)-c(j,i)
1555 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1559 if (itype(i).ne.10) then
1561 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1562 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1570 c----------------------------------------------------------------------------
1571 subroutine read_threadbase
1572 implicit real*8 (a-h,o-z)
1573 include 'DIMENSIONS'
1574 include 'COMMON.IOUNITS'
1575 include 'COMMON.GEO'
1576 include 'COMMON.VAR'
1577 include 'COMMON.INTERACT'
1578 include 'COMMON.LOCAL'
1579 include 'COMMON.NAMES'
1580 include 'COMMON.CHAIN'
1581 include 'COMMON.FFIELD'
1582 include 'COMMON.SBRIDGE'
1583 include 'COMMON.HEADER'
1584 include 'COMMON.CONTROL'
1585 include 'COMMON.DBASE'
1586 include 'COMMON.THREAD'
1587 include 'COMMON.TIME1'
1588 C Read pattern database for threading.
1589 read (icbase,*) nseq
1591 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1592 & nres_base(2,i),nres_base(3,i)
1593 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1595 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1596 c & nres_base(2,i),nres_base(3,i)
1597 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1601 if (weidis.eq.0.0D0) weidis=0.1D0
1610 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1611 write (iout,'(a,i5)') 'nexcl: ',nexcl
1612 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1615 c------------------------------------------------------------------------------
1616 subroutine setup_var
1617 implicit real*8 (a-h,o-z)
1618 include 'DIMENSIONS'
1619 include 'COMMON.IOUNITS'
1620 include 'COMMON.GEO'
1621 include 'COMMON.VAR'
1622 include 'COMMON.INTERACT'
1623 include 'COMMON.LOCAL'
1624 include 'COMMON.NAMES'
1625 include 'COMMON.CHAIN'
1626 include 'COMMON.FFIELD'
1627 include 'COMMON.SBRIDGE'
1628 include 'COMMON.HEADER'
1629 include 'COMMON.CONTROL'
1630 include 'COMMON.DBASE'
1631 include 'COMMON.THREAD'
1632 include 'COMMON.TIME1'
1633 C Set up variable list.
1639 if (itype(i).ne.10) then
1641 ialph(i,1)=nvar+nside
1645 if (indphi.gt.0) then
1647 else if (indback.gt.0) then
1652 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1655 c----------------------------------------------------------------------------
1656 subroutine gen_dist_constr
1657 C Generate CA distance constraints.
1658 implicit real*8 (a-h,o-z)
1659 include 'DIMENSIONS'
1660 include 'COMMON.IOUNITS'
1661 include 'COMMON.GEO'
1662 include 'COMMON.VAR'
1663 include 'COMMON.INTERACT'
1664 include 'COMMON.LOCAL'
1665 include 'COMMON.NAMES'
1666 include 'COMMON.CHAIN'
1667 include 'COMMON.FFIELD'
1668 include 'COMMON.SBRIDGE'
1669 include 'COMMON.HEADER'
1670 include 'COMMON.CONTROL'
1671 include 'COMMON.DBASE'
1672 include 'COMMON.THREAD'
1673 include 'COMMON.TIME1'
1674 dimension itype_pdb(maxres)
1675 common /pizda/ itype_pdb
1677 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1678 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1679 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1681 do i=nstart_sup,nstart_sup+nsup-1
1682 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1683 cd & ' seq_pdb', restyp(itype_pdb(i))
1684 do j=i+2,nstart_sup+nsup-1
1686 ihpb(nhpb)=i+nstart_seq-nstart_sup
1687 jhpb(nhpb)=j+nstart_seq-nstart_sup
1689 dhpb(nhpb)=dist(i,j)
1692 cd write (iout,'(a)') 'Distance constraints:'
1697 cd if (ii.gt.nres) then
1702 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1703 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1704 cd & dhpb(i),forcon(i)
1708 c----------------------------------------------------------------------------
1710 implicit real*8 (a-h,o-z)
1711 include 'DIMENSIONS'
1712 include 'COMMON.MAP'
1713 include 'COMMON.IOUNITS'
1714 character*3 angid(4) /'THE','PHI','ALP','OME'/
1715 character*80 mapcard,ucase
1717 read (inp,'(a)') mapcard
1718 mapcard=ucase(mapcard)
1719 if (index(mapcard,'PHI').gt.0) then
1721 else if (index(mapcard,'THE').gt.0) then
1723 else if (index(mapcard,'ALP').gt.0) then
1725 else if (index(mapcard,'OME').gt.0) then
1728 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1729 stop 'Error - illegal variable spec in MAP card.'
1731 call readi (mapcard,'RES1',res1(imap),0)
1732 call readi (mapcard,'RES2',res2(imap),0)
1733 if (res1(imap).eq.0) then
1734 res1(imap)=res2(imap)
1735 else if (res2(imap).eq.0) then
1736 res2(imap)=res1(imap)
1738 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1740 & 'Error - illegal definition of variable group in MAP.'
1741 stop 'Error - illegal definition of variable group in MAP.'
1743 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1744 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1745 call readi(mapcard,'NSTEP',nstep(imap),0)
1746 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1748 & 'Illegal boundary and/or step size specification in MAP.'
1749 stop 'Illegal boundary and/or step size specification in MAP.'
1754 c----------------------------------------------------------------------------
1755 csa subroutine csaread
1756 csa implicit real*8 (a-h,o-z)
1757 csa include 'DIMENSIONS'
1758 csa include 'COMMON.IOUNITS'
1759 csa include 'COMMON.GEO'
1760 csa include 'COMMON.CSA'
1761 csa include 'COMMON.BANK'
1762 csa include 'COMMON.CONTROL'
1763 csa character*80 ucase
1764 csa character*620 mcmcard
1765 csa call card_concat(mcmcard)
1767 csa call readi(mcmcard,'NCONF',nconf,50)
1768 csa call readi(mcmcard,'NADD',nadd,0)
1769 csa call readi(mcmcard,'JSTART',jstart,1)
1770 csa call readi(mcmcard,'JEND',jend,1)
1771 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1772 csa call readi(mcmcard,'N0',n0,1)
1773 csa call readi(mcmcard,'N1',n1,6)
1774 csa call readi(mcmcard,'N2',n2,4)
1775 csa call readi(mcmcard,'N3',n3,0)
1776 csa call readi(mcmcard,'N4',n4,0)
1777 csa call readi(mcmcard,'N5',n5,0)
1778 csa call readi(mcmcard,'N6',n6,10)
1779 csa call readi(mcmcard,'N7',n7,0)
1780 csa call readi(mcmcard,'N8',n8,0)
1781 csa call readi(mcmcard,'N9',n9,0)
1782 csa call readi(mcmcard,'N14',n14,0)
1783 csa call readi(mcmcard,'N15',n15,0)
1784 csa call readi(mcmcard,'N16',n16,0)
1785 csa call readi(mcmcard,'N17',n17,0)
1786 csa call readi(mcmcard,'N18',n18,0)
1788 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1790 csa call readi(mcmcard,'NDIFF',ndiff,2)
1791 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1792 csa call readi(mcmcard,'IS1',is1,1)
1793 csa call readi(mcmcard,'IS2',is2,8)
1794 csa call readi(mcmcard,'NRAN0',nran0,4)
1795 csa call readi(mcmcard,'NRAN1',nran1,2)
1796 csa call readi(mcmcard,'IRR',irr,1)
1797 csa call readi(mcmcard,'NSEED',nseed,20)
1798 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1799 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1800 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1801 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1802 csa call readi(mcmcard,'ICMAX',icmax,3)
1803 csa call readi(mcmcard,'IRESTART',irestart,0)
1804 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1807 csa call reada(mcmcard,'DELE',dele,20.0d0)
1808 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1809 csa call readi(mcmcard,'IREF',iref,0)
1810 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1811 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1812 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1813 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1814 csa write (iout,*) "NCONF_IN",nconf_in
1817 c----------------------------------------------------------------------------
1818 cfmc subroutine mcmfread
1819 cfmc implicit real*8 (a-h,o-z)
1820 cfmc include 'DIMENSIONS'
1821 cfmc include 'COMMON.MCMF'
1822 cfmc include 'COMMON.IOUNITS'
1823 cfmc include 'COMMON.GEO'
1824 cfmc character*80 ucase
1825 cfmc character*620 mcmcard
1826 cfmc call card_concat(mcmcard)
1828 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1829 cfmc write(iout,*)'MAXRANT=',maxrant
1830 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1831 cfmc write(iout,*)'MAXFAM=',maxfam
1832 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1833 cfmc write(iout,*)'NNET1=',nnet1
1834 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1835 cfmc write(iout,*)'NNET2=',nnet2
1836 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1837 cfmc write(iout,*)'NNET3=',nnet3
1838 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1839 cfmc write(iout,*)'ILASTT=',ilastt
1840 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1841 cfmc write(iout,*)'MAXSTR=',maxstr
1842 cfmc maxstr_f=maxstr/maxfam
1843 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1844 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1845 cfmc write(iout,*)'NMCMF=',nmcmf
1846 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1847 cfmc write(iout,*)'IFOCUS=',ifocus
1848 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1849 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1850 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1851 cfmc write(iout,*)'INTPRT=',intprt
1852 cfmc call readi(mcmcard,'IPRT',iprt,100)
1853 cfmc write(iout,*)'IPRT=',iprt
1854 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1855 cfmc write(iout,*)'IMAXTR=',imaxtr
1856 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1857 cfmc write(iout,*)'MAXEVEN=',maxeven
1858 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1859 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1860 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1861 cfmc write(iout,*)'INIMIN=',inimin
1862 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1863 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1864 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1865 cfmc write(iout,*)'NTHREAD=',nthread
1866 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1867 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1868 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1869 cfmc write(iout,*)'MAXPERT=',maxpert
1870 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1871 cfmc write(iout,*)'IRMSD=',irmsd
1872 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1873 cfmc write(iout,*)'DENEMIN=',denemin
1874 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1875 cfmc write(iout,*)'RCUT1S=',rcut1s
1876 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1877 cfmc write(iout,*)'RCUT1E=',rcut1e
1878 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1879 cfmc write(iout,*)'RCUT2S=',rcut2s
1880 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1881 cfmc write(iout,*)'RCUT2E=',rcut2e
1882 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1883 cfmc write(iout,*)'DPERT1=',d_pert1
1884 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1885 cfmc write(iout,*)'DPERT1A=',d_pert1a
1886 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1887 cfmc write(iout,*)'DPERT2=',d_pert2
1888 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1889 cfmc write(iout,*)'DPERT2A=',d_pert2a
1890 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1891 cfmc write(iout,*)'DPERT2B=',d_pert2b
1892 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1893 cfmc write(iout,*)'DPERT2C=',d_pert2c
1894 cfmc d_pert1=deg2rad*d_pert1
1895 cfmc d_pert1a=deg2rad*d_pert1a
1896 cfmc d_pert2=deg2rad*d_pert2
1897 cfmc d_pert2a=deg2rad*d_pert2a
1898 cfmc d_pert2b=deg2rad*d_pert2b
1899 cfmc d_pert2c=deg2rad*d_pert2c
1900 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1901 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1902 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1903 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1904 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1905 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1906 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1907 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1908 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1909 cfmc write(iout,*)'RCUTINI=',rcutini
1910 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1911 cfmc write(iout,*)'GRAT=',grat
1912 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1913 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1917 c----------------------------------------------------------------------------
1919 implicit real*8 (a-h,o-z)
1920 include 'DIMENSIONS'
1921 include 'COMMON.MCM'
1922 include 'COMMON.MCE'
1923 include 'COMMON.IOUNITS'
1925 character*320 mcmcard
1926 call card_concat(mcmcard)
1927 call readi(mcmcard,'MAXACC',maxacc,100)
1928 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1929 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1930 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1931 call readi(mcmcard,'MAXREPM',maxrepm,200)
1932 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1933 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1934 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1935 call reada(mcmcard,'E_UP',e_up,5.0D0)
1936 call reada(mcmcard,'DELTE',delte,0.1D0)
1937 call readi(mcmcard,'NSWEEP',nsweep,5)
1938 call readi(mcmcard,'NSTEPH',nsteph,0)
1939 call readi(mcmcard,'NSTEPC',nstepc,0)
1940 call reada(mcmcard,'TMIN',tmin,298.0D0)
1941 call reada(mcmcard,'TMAX',tmax,298.0D0)
1942 call readi(mcmcard,'NWINDOW',nwindow,0)
1943 call readi(mcmcard,'PRINT_MC',print_mc,0)
1944 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1945 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1946 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1947 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1948 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1949 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1950 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1951 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1952 if (nwindow.gt.0) then
1953 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1955 winlen(i)=winend(i)-winstart(i)+1
1958 if (tmax.lt.tmin) tmax=tmin
1959 if (tmax.eq.tmin) then
1963 if (nstepc.gt.0 .and. nsteph.gt.0) then
1964 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1965 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1967 C Probabilities of different move types
1968 sumpro_type(0)=0.0D0
1969 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1970 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1971 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1972 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1973 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1974 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1975 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1977 print *,'i',i,' sumprotype',sumpro_type(i)
1978 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1979 print *,'i',i,' sumprotype',sumpro_type(i)
1983 c----------------------------------------------------------------------------
1984 subroutine read_minim
1985 implicit real*8 (a-h,o-z)
1986 include 'DIMENSIONS'
1987 include 'COMMON.MINIM'
1988 include 'COMMON.IOUNITS'
1990 character*320 minimcard
1991 call card_concat(minimcard)
1992 call readi(minimcard,'MAXMIN',maxmin,2000)
1993 call readi(minimcard,'MAXFUN',maxfun,5000)
1994 call readi(minimcard,'MINMIN',minmin,maxmin)
1995 call readi(minimcard,'MINFUN',minfun,maxmin)
1996 call reada(minimcard,'TOLF',tolf,1.0D-2)
1997 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1998 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1999 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
2000 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
2001 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2002 & 'Options in energy minimization:'
2003 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
2004 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
2005 & 'MinMin:',MinMin,' MinFun:',MinFun,
2006 & ' TolF:',TolF,' RTolF:',RTolF
2009 c----------------------------------------------------------------------------
2010 subroutine read_angles(kanal,*)
2011 implicit real*8 (a-h,o-z)
2012 include 'DIMENSIONS'
2013 include 'COMMON.GEO'
2014 include 'COMMON.VAR'
2015 include 'COMMON.CHAIN'
2016 include 'COMMON.IOUNITS'
2017 include 'COMMON.CONTROL'
2018 c Read angles from input
2020 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
2021 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
2022 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
2023 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
2026 c 9/7/01 avoid 180 deg valence angle
2027 if (theta(i).gt.179.99d0) theta(i)=179.99d0
2029 theta(i)=deg2rad*theta(i)
2030 phi(i)=deg2rad*phi(i)
2031 alph(i)=deg2rad*alph(i)
2032 omeg(i)=deg2rad*omeg(i)
2037 c----------------------------------------------------------------------------
2038 subroutine reada(rekord,lancuch,wartosc,default)
2040 character*(*) rekord,lancuch
2041 double precision wartosc,default
2044 iread=index(rekord,lancuch)
2045 if (iread.eq.0) then
2049 iread=iread+ilen(lancuch)+1
2050 read (rekord(iread:),*,err=10,end=10) wartosc
2055 c----------------------------------------------------------------------------
2056 subroutine readi(rekord,lancuch,wartosc,default)
2058 character*(*) rekord,lancuch
2059 integer wartosc,default
2062 iread=index(rekord,lancuch)
2063 if (iread.eq.0) then
2067 iread=iread+ilen(lancuch)+1
2068 read (rekord(iread:),*,err=10,end=10) wartosc
2073 c----------------------------------------------------------------------------
2074 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2077 integer tablica(dim),default
2078 character*(*) rekord,lancuch
2085 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2086 if (iread.eq.0) return
2087 iread=iread+ilen(lancuch)+1
2088 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2091 c----------------------------------------------------------------------------
2092 subroutine multreada(rekord,lancuch,tablica,dim,default)
2095 double precision tablica(dim),default
2096 character*(*) rekord,lancuch
2103 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2104 if (iread.eq.0) return
2105 iread=iread+ilen(lancuch)+1
2106 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2109 c----------------------------------------------------------------------------
2110 subroutine openunits
2111 implicit real*8 (a-h,o-z)
2112 include 'DIMENSIONS'
2115 character*16 form,nodename
2118 include 'COMMON.SETUP'
2119 include 'COMMON.IOUNITS'
2121 include 'COMMON.CONTROL'
2122 integer lenpre,lenpot,ilen,lentmp
2124 character*3 out1file_text,ucase
2127 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2128 call getenv_loc("PREFIX",prefix)
2130 call getenv_loc("POT",pot)
2131 call getenv_loc("DIRTMP",tmpdir)
2132 call getenv_loc("CURDIR",curdir)
2133 call getenv_loc("OUT1FILE",out1file_text)
2134 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2135 out1file_text=ucase(out1file_text)
2136 if (out1file_text(1:1).eq."Y") then
2139 out1file=fg_rank.gt.0
2144 if (lentmp.gt.0) then
2145 write (*,'(80(1h!))')
2146 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2147 write (*,'(80(1h!))')
2148 write (*,*)"All output files will be on node /tmp directory."
2150 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2151 if (me.eq.king) then
2152 write (*,*) "The master node is ",nodename
2153 else if (fg_rank.eq.0) then
2154 write (*,*) "I am the CG slave node ",nodename
2156 write (*,*) "I am the FG slave node ",nodename
2159 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2160 lenpre = lentmp+lenpre+1
2162 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2163 C Get the names and open the input files
2164 #if defined(WINIFL) || defined(WINPGI)
2165 open(1,file=pref_orig(:ilen(pref_orig))//
2166 & '.inp',status='old',readonly,shared)
2167 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2168 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2169 C Get parameter filenames and open the parameter files.
2170 call getenv_loc('BONDPAR',bondname)
2171 open (ibond,file=bondname,status='old',readonly,shared)
2172 call getenv_loc('THETPAR',thetname)
2173 open (ithep,file=thetname,status='old',readonly,shared)
2175 call getenv_loc('THETPARPDB',thetname_pdb)
2176 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2178 call getenv_loc('ROTPAR',rotname)
2179 open (irotam,file=rotname,status='old',readonly,shared)
2181 call getenv_loc('ROTPARPDB',rotname_pdb)
2182 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2184 call getenv_loc('TORPAR',torname)
2185 open (itorp,file=torname,status='old',readonly,shared)
2186 call getenv_loc('TORDPAR',tordname)
2187 open (itordp,file=tordname,status='old',readonly,shared)
2188 call getenv_loc('FOURIER',fouriername)
2189 open (ifourier,file=fouriername,status='old',readonly,shared)
2190 call getenv_loc('ELEPAR',elename)
2191 open (ielep,file=elename,status='old',readonly,shared)
2192 call getenv_loc('SIDEPAR',sidename)
2193 open (isidep,file=sidename,status='old',readonly,shared)
2194 #elif (defined CRAY) || (defined AIX)
2195 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2197 c print *,"Processor",myrank," opened file 1"
2198 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2199 c print *,"Processor",myrank," opened file 9"
2200 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2201 C Get parameter filenames and open the parameter files.
2202 call getenv_loc('BONDPAR',bondname)
2203 open (ibond,file=bondname,status='old',action='read')
2204 c print *,"Processor",myrank," opened file IBOND"
2205 call getenv_loc('THETPAR',thetname)
2206 open (ithep,file=thetname,status='old',action='read')
2207 c print *,"Processor",myrank," opened file ITHEP"
2209 call getenv_loc('THETPARPDB',thetname_pdb)
2210 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2212 call getenv_loc('ROTPAR',rotname)
2213 open (irotam,file=rotname,status='old',action='read')
2214 c print *,"Processor",myrank," opened file IROTAM"
2216 call getenv_loc('ROTPARPDB',rotname_pdb)
2217 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2219 call getenv_loc('TORPAR',torname)
2220 open (itorp,file=torname,status='old',action='read')
2221 c print *,"Processor",myrank," opened file ITORP"
2222 call getenv_loc('TORDPAR',tordname)
2223 open (itordp,file=tordname,status='old',action='read')
2224 c print *,"Processor",myrank," opened file ITORDP"
2225 call getenv_loc('SCCORPAR',sccorname)
2226 open (isccor,file=sccorname,status='old',action='read')
2227 c print *,"Processor",myrank," opened file ISCCOR"
2228 call getenv_loc('FOURIER',fouriername)
2229 open (ifourier,file=fouriername,status='old',action='read')
2230 c print *,"Processor",myrank," opened file IFOURIER"
2231 call getenv_loc('ELEPAR',elename)
2232 open (ielep,file=elename,status='old',action='read')
2233 c print *,"Processor",myrank," opened file IELEP"
2234 call getenv_loc('SIDEPAR',sidename)
2235 open (isidep,file=sidename,status='old',action='read')
2236 c print *,"Processor",myrank," opened file ISIDEP"
2237 c print *,"Processor",myrank," opened parameter files"
2239 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2240 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2241 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2242 C Get parameter filenames and open the parameter files.
2243 call getenv_loc('BONDPAR',bondname)
2244 open (ibond,file=bondname,status='old')
2245 call getenv_loc('THETPAR',thetname)
2246 open (ithep,file=thetname,status='old')
2248 call getenv_loc('THETPARPDB',thetname_pdb)
2249 open (ithep_pdb,file=thetname_pdb,status='old')
2251 call getenv_loc('ROTPAR',rotname)
2252 open (irotam,file=rotname,status='old')
2254 call getenv_loc('ROTPARPDB',rotname_pdb)
2255 open (irotam_pdb,file=rotname_pdb,status='old')
2257 call getenv_loc('TORPAR',torname)
2258 open (itorp,file=torname,status='old')
2259 call getenv_loc('TORDPAR',tordname)
2260 open (itordp,file=tordname,status='old')
2261 call getenv_loc('SCCORPAR',sccorname)
2262 open (isccor,file=sccorname,status='old')
2263 call getenv_loc('FOURIER',fouriername)
2264 open (ifourier,file=fouriername,status='old')
2265 call getenv_loc('ELEPAR',elename)
2266 open (ielep,file=elename,status='old')
2267 call getenv_loc('SIDEPAR',sidename)
2268 open (isidep,file=sidename,status='old')
2270 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2272 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2273 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2274 C Get parameter filenames and open the parameter files.
2275 call getenv_loc('BONDPAR',bondname)
2276 open (ibond,file=bondname,status='old',action='read')
2277 call getenv_loc('THETPAR',thetname)
2278 open (ithep,file=thetname,status='old',action='read')
2280 call getenv_loc('THETPARPDB',thetname_pdb)
2281 print *,"thetname_pdb ",thetname_pdb
2282 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2283 print *,ithep_pdb," opened"
2285 call getenv_loc('ROTPAR',rotname)
2286 open (irotam,file=rotname,status='old',action='read')
2288 call getenv_loc('ROTPARPDB',rotname_pdb)
2289 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2291 call getenv_loc('TORPAR',torname)
2292 open (itorp,file=torname,status='old',action='read')
2293 call getenv_loc('TORDPAR',tordname)
2294 open (itordp,file=tordname,status='old',action='read')
2295 call getenv_loc('SCCORPAR',sccorname)
2296 open (isccor,file=sccorname,status='old',action='read')
2297 call getenv_loc('FOURIER',fouriername)
2298 open (ifourier,file=fouriername,status='old',action='read')
2299 call getenv_loc('ELEPAR',elename)
2300 open (ielep,file=elename,status='old',action='read')
2301 call getenv_loc('SIDEPAR',sidename)
2302 open (isidep,file=sidename,status='old',action='read')
2306 C 8/9/01 In the newest version SCp interaction constants are read from a file
2307 C Use -DOLDSCP to use hard-coded constants instead.
2309 call getenv_loc('SCPPAR',scpname)
2310 #if defined(WINIFL) || defined(WINPGI)
2311 open (iscpp,file=scpname,status='old',readonly,shared)
2312 #elif (defined CRAY) || (defined AIX)
2313 open (iscpp,file=scpname,status='old',action='read')
2315 open (iscpp,file=scpname,status='old')
2317 open (iscpp,file=scpname,status='old',action='read')
2320 call getenv_loc('PATTERN',patname)
2321 #if defined(WINIFL) || defined(WINPGI)
2322 open (icbase,file=patname,status='old',readonly,shared)
2323 #elif (defined CRAY) || (defined AIX)
2324 open (icbase,file=patname,status='old',action='read')
2326 open (icbase,file=patname,status='old')
2328 open (icbase,file=patname,status='old',action='read')
2331 C Open output file only for CG processes
2332 c print *,"Processor",myrank," fg_rank",fg_rank
2333 if (fg_rank.eq.0) then
2335 if (nodes.eq.1) then
2338 npos = dlog10(dfloat(nodes-1))+1
2340 if (npos.lt.3) npos=3
2341 write (liczba,'(i1)') npos
2342 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2344 write (liczba,form) me
2345 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2346 & liczba(:ilen(liczba))
2347 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2349 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2351 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2352 & liczba(:ilen(liczba))//'.mol2'
2353 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//
2354 & liczba(:ilen(liczba))//'.x'
2355 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2356 & liczba(:ilen(liczba))//'.stat'
2358 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2359 & //liczba(:ilen(liczba))//'.stat')
2360 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2363 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2364 & liczba(:ilen(liczba))//'.const'
2369 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2370 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2371 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2372 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2373 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x'
2374 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2376 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2378 rest2name=prefix(:ilen(prefix))//'.rst'
2380 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2383 #if defined(AIX) || defined(PGI)
2384 if (me.eq.king .or. .not. out1file)
2385 & open(iout,file=outname,status='unknown')
2388 if (fg_rank.gt.0) then
2389 write (liczba,'(i3.3)') myrank/nfgtasks
2390 write (ll,'(bz,i3.3)') fg_rank
2391 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2397 open(igeom,file=intname,status='unknown',position='append')
2398 open(ipdb,file=pdbname,status='unknown')
2399 open(imol2,file=mol2name,status='unknown')
2400 open(istat,file=statname,status='unknown',position='append')
2402 c1out open(iout,file=outname,status='unknown')
2405 if (me.eq.king .or. .not.out1file)
2406 & open(iout,file=outname,status='unknown')
2409 if (fg_rank.gt.0) then
2410 print "Processor",fg_rank," opening output file"
2411 write (liczba,'(i3.3)') myrank/nfgtasks
2412 write (ll,'(bz,i3.3)') fg_rank
2413 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2419 open(igeom,file=intname,status='unknown',access='append')
2420 open(ipdb,file=pdbname,status='unknown')
2421 open(imol2,file=mol2name,status='unknown')
2422 open(istat,file=statname,status='unknown',access='append')
2424 c1out open(iout,file=outname,status='unknown')
2427 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2428 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2429 csa csa_history=prefix(:lenpre)//'.CSA.history'
2430 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2431 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2432 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2433 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2434 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2435 csa csa_int=prefix(:lenpre)//'.int'
2436 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2437 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2438 csa csa_in=prefix(:lenpre)//'.CSA.in'
2439 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2442 write (iout,'(80(1h-))')
2443 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2444 write (iout,'(80(1h-))')
2445 write (iout,*) "Input file : ",
2446 & pref_orig(:ilen(pref_orig))//'.inp'
2447 write (iout,*) "Output file : ",
2448 & outname(:ilen(outname))
2450 write (iout,*) "Sidechain potential file : ",
2451 & sidename(:ilen(sidename))
2453 write (iout,*) "SCp potential file : ",
2454 & scpname(:ilen(scpname))
2456 write (iout,*) "Electrostatic potential file : ",
2457 & elename(:ilen(elename))
2458 write (iout,*) "Cumulant coefficient file : ",
2459 & fouriername(:ilen(fouriername))
2460 write (iout,*) "Torsional parameter file : ",
2461 & torname(:ilen(torname))
2462 write (iout,*) "Double torsional parameter file : ",
2463 & tordname(:ilen(tordname))
2464 write (iout,*) "SCCOR parameter file : ",
2465 & sccorname(:ilen(sccorname))
2466 write (iout,*) "Bond & inertia constant file : ",
2467 & bondname(:ilen(bondname))
2468 write (iout,*) "Bending parameter file : ",
2469 & thetname(:ilen(thetname))
2470 write (iout,*) "Rotamer parameter file : ",
2471 & rotname(:ilen(rotname))
2472 write (iout,*) "Threading database : ",
2473 & patname(:ilen(patname))
2475 &write (iout,*)" DIRTMP : ",
2477 write (iout,'(80(1h-))')
2481 c----------------------------------------------------------------------------
2482 subroutine card_concat(card)
2483 implicit real*8 (a-h,o-z)
2484 include 'DIMENSIONS'
2485 include 'COMMON.IOUNITS'
2487 character*80 karta,ucase
2489 read (inp,'(a)') karta
2492 do while (karta(80:80).eq.'&')
2493 card=card(:ilen(card)+1)//karta(:79)
2494 read (inp,'(a)') karta
2497 card=card(:ilen(card)+1)//karta
2500 c----------------------------------------------------------------------------------
2502 implicit real*8 (a-h,o-z)
2503 include 'DIMENSIONS'
2504 include 'COMMON.CHAIN'
2505 include 'COMMON.IOUNITS'
2507 include 'COMMON.CONTROL'
2508 open(irest2,file=rest2name,status='unknown')
2509 read(irest2,*) totT,EK,potE,totE,t_bath
2511 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2514 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2516 if(usampl.or.homol_nset.gt.1) then
2517 read (irest2,*) iset
2522 c---------------------------------------------------------------------------------
2523 subroutine read_fragments
2524 implicit real*8 (a-h,o-z)
2525 include 'DIMENSIONS'
2529 include 'COMMON.SETUP'
2530 include 'COMMON.CHAIN'
2531 include 'COMMON.IOUNITS'
2533 include 'COMMON.CONTROL'
2535 read(inp,*) nset,nfrag,npair,nfrag_back
2536 if(me.eq.king.or..not.out1file)
2537 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2538 & " nfrag_back",nfrag_back
2540 read(inp,*) mset(iset1)
2542 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2544 if(me.eq.king.or..not.out1file)
2545 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2546 & ifrag(2,i,iset1), qinfrag(i,iset1)
2549 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2551 if(me.eq.king.or..not.out1file)
2552 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2553 & ipair(2,i,iset1), qinpair(i,iset1)
2556 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2557 & wfrag_back(3,i,iset1),
2558 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2559 if(me.eq.king.or..not.out1file)
2560 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2561 & wfrag_back(2,i,iset1),
2562 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2563 & ifrag_back(2,i,iset1)
2568 c-------------------------------------------------------------------------------
2569 subroutine read_dist_constr
2570 implicit real*8 (a-h,o-z)
2571 include 'DIMENSIONS'
2575 include 'COMMON.SETUP'
2576 include 'COMMON.CONTROL'
2577 include 'COMMON.CHAIN'
2578 include 'COMMON.IOUNITS'
2579 include 'COMMON.SBRIDGE'
2580 integer ifrag_(2,100),ipair_(2,100)
2581 double precision wfrag_(100),wpair_(100)
2582 character*500 controlcard
2583 logical normalize,next
2585 double precision xlink(4,0:4) /
2587 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2588 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2589 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2590 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2591 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2592 c print *, "WCHODZE"
2593 write (iout,*) "Calling read_dist_constr"
2594 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2600 call card_concat(controlcard)
2601 next = index(controlcard,"NEXT").gt.0
2602 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2603 write (iout,*) "restr_type",restr_type
2604 call readi(controlcard,"NFRAG",nfrag_,0)
2605 call readi(controlcard,"NFRAG",nfrag_,0)
2606 call readi(controlcard,"NPAIR",npair_,0)
2607 call readi(controlcard,"NDIST",ndist_,0)
2608 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2609 if (restr_type.eq.10)
2610 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2611 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2612 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2613 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2614 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2615 normalize = index(controlcard,"NORMALIZE").gt.0
2616 write (iout,*) "WBOLTZD",wboltzd
2617 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2618 c write (iout,*) "IFRAG"
2620 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2622 c write (iout,*) "IPAIR"
2624 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2626 if (nfrag_.gt.0) write (iout,*)
2627 & "Distance restraints as generated from reference structure"
2629 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2630 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2631 & ifrag_(2,i)=nstart_sup+nsup-1
2632 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2634 if (wfrag_(i).eq.0.0d0) cycle
2635 do j=ifrag_(1,i),ifrag_(2,i)-1
2636 do k=j+1,ifrag_(2,i)
2637 c write (iout,*) "j",j," k",k
2639 if (restr_type.eq.1) then
2645 forcon(nhpb)=wfrag_(i)
2646 else if (constr_dist.eq.2) then
2647 if (ddjk.le.dist_cut) then
2653 forcon(nhpb)=wfrag_(i)
2655 else if (restr_type.eq.3) then
2661 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2664 if (.not.out1file .or. me.eq.king)
2665 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2666 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2668 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2669 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2675 if (wpair_(i).eq.0.0d0) cycle
2683 do j=ifrag_(1,ii),ifrag_(2,ii)
2684 do k=ifrag_(1,jj),ifrag_(2,jj)
2685 if (restr_type.eq.1) then
2691 forcon(nhpb)=wfrag_(i)
2692 else if (constr_dist.eq.2) then
2693 if (ddjk.le.dist_cut) then
2699 forcon(nhpb)=wfrag_(i)
2701 else if (restr_type.eq.3) then
2707 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2710 if (.not.out1file .or. me.eq.king)
2711 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
2712 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2714 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
2715 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2722 write (iout,*) "Distance restraints as read from input"
2724 if (restr_type.eq.11) then
2725 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2726 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2727 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2728 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2730 irestr_type(nhpb)=11
2732 if (.not.out1file .or. me.eq.king)
2733 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2734 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2735 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2737 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2738 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2739 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2741 if (ibecarb(nhpb).gt.0) then
2742 ihpb(nhpb)=ihpb(nhpb)+nres
2743 jhpb(nhpb)=jhpb(nhpb)+nres
2745 else if (constr_dist.eq.10) then
2746 c Cross-lonk Markov-like potential
2747 call card_concat(controlcard)
2748 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2749 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2751 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2752 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2753 if (index(controlcard,"ZL").gt.0) then
2755 else if (index(controlcard,"ADH").gt.0) then
2757 else if (index(controlcard,"PDH").gt.0) then
2759 else if (index(controlcard,"DSS").gt.0) then
2764 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2765 & xlink(1,link_type))
2766 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2767 & xlink(2,link_type))
2768 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2769 & xlink(3,link_type))
2770 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2771 & xlink(4,link_type))
2772 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2773 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2774 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2775 if (forcon(nhpb+1).le.0.0d0 .or.
2776 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2778 irestr_type(nhpb)=10
2779 if (ibecarb(nhpb).gt.0) then
2780 ihpb(nhpb)=ihpb(nhpb)+nres
2781 jhpb(nhpb)=jhpb(nhpb)+nres
2784 if (.not.out1file .or. me.eq.king)
2785 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2786 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2787 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2790 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2791 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2792 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2797 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2798 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2799 if (forcon(nhpb+1).gt.0.0d0) then
2801 if (dhpb1(nhpb).eq.0.0d0) then
2806 if (ibecarb(nhpb).gt.0) then
2807 ihpb(nhpb)=ihpb(nhpb)+nres
2808 jhpb(nhpb)=jhpb(nhpb)+nres
2810 if (dhpb(nhpb).eq.0.0d0)
2811 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2814 if (.not.out1file .or. me.eq.king)
2815 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2816 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2818 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2819 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2822 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2823 C if (forcon(nhpb+1).gt.0.0d0) then
2825 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2833 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2834 & fordepthmax=fordepth(i)
2837 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2842 c-------------------------------------------------------------------------------
2844 subroutine read_constr_homology
2846 include 'DIMENSIONS'
2850 include 'COMMON.SETUP'
2851 include 'COMMON.CONTROL'
2852 include 'COMMON.CHAIN'
2853 include 'COMMON.IOUNITS'
2855 include 'COMMON.GEO'
2856 include 'COMMON.INTERACT'
2857 include 'COMMON.NAMES'
2859 c For new homol impl
2861 include 'COMMON.VAR'
2864 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2866 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2867 c & sigma_odl_temp(maxres,maxres,max_template)
2869 character*24 model_ki_dist, model_ki_angle
2870 character*500 controlcard
2871 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2872 logical lprn /.true./
2877 c FP - Nov. 2014 Temporary specifications for new vars
2879 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
2880 double precision, dimension (max_template,maxres) :: rescore
2881 double precision, dimension (max_template,maxres) :: rescore2
2882 character*24 pdbfile,tpl_k_rescore
2883 c -----------------------------------------------------------------
2884 c Reading multiple PDB ref structures and calculation of retraints
2885 c not using pre-computed ones stored in files model_ki_{dist,angle}
2887 c -----------------------------------------------------------------
2890 c Alternative: reading from input
2891 call card_concat(controlcard)
2892 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2893 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2894 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2895 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2896 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2897 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2898 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2899 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2900 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2901 if(.not.read2sigma.and.start_from_model) then
2902 write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2903 start_from_model=.false.
2905 if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON'
2906 if(start_from_model .and. rest) then
2907 write(iout,*) 'START_FROM_MODELS is OFF'
2908 write(iout,*) 'remove restart keyword from input'
2910 if (homol_nset.gt.1)then
2911 call card_concat(controlcard)
2912 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2913 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2914 write(iout,*) "iset homology_weight "
2916 write(iout,*) i,waga_homology(i)
2919 iset=mod(kolor,homol_nset)+1
2922 waga_homology(1)=1.0
2925 cd write (iout,*) "nnt",nnt," nct",nct
2932 write(iout,*) 'nnt=',nnt,'nct=',nct
2935 do k=1,constr_homology
2948 if (read_homol_frag) then
2949 call read_klapaucjusz
2952 do k=1,constr_homology
2954 read(inp,'(a)') pdbfile
2955 c Next stament causes error upon compilation (?)
2956 c if(me.eq.king.or. .not. out1file)
2957 c write (iout,'(2a)') 'PDB data will be read from file ',
2958 c & pdbfile(:ilen(pdbfile))
2959 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2960 & pdbfile(:ilen(pdbfile))
2961 open(ipdbin,file=pdbfile,status='old',err=33)
2963 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2964 & pdbfile(:ilen(pdbfile))
2967 c print *,'Begin reading pdb data'
2969 c Files containing res sim or local scores (former containing sigmas)
2972 write(kic2,'(bz,i2.2)') k
2974 tpl_k_rescore="template"//kic2//".sco"
2977 if (read2sigma) then
2978 call readpdb_template(k)
2983 c Distance restraints
2986 C Copy the coordinates from reference coordinates (?)
2990 c write (iout,*) "c(",j,i,") =",c(j,i)
2994 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2996 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2997 open (ientin,file=tpl_k_rescore,status='old')
2998 if (nnt.gt.1) rescore(k,1)=0.0d0
2999 do irec=nnt,nct ! loop for reading res sim
3000 if (read2sigma) then
3001 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3004 idomain(k,i_tmp)=idomain_tmp
3005 rescore(k,i_tmp)=rescore_tmp
3006 rescore2(k,i_tmp)=rescore2_tmp
3007 write(iout,'(a7,i5,2f10.5,i5)') "rescore",
3008 & i_tmp,rescore2_tmp,rescore_tmp,
3012 read (ientin,*,end=1401) rescore_tmp
3014 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3015 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3016 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3021 if (waga_dist.ne.0.0d0) then
3029 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3030 c write (iout,*) k,i,j,distal,dist2_cut
3032 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3033 & .and. distal.le.dist2_cut ) then
3039 c write (iout,*) "k",k
3040 c write (iout,*) "i",i," j",j," constr_homology",
3045 if (read2sigma) then
3048 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3050 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3051 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3052 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3054 if (odl(k,ii).le.dist_cut) then
3055 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3058 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3059 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3061 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3062 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3066 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3069 l_homo(k,ii)=.false.
3076 c Theta, dihedral and SC retraints
3078 if (waga_angle.gt.0.0d0) then
3079 c open (ientin,file=tpl_k_sigma_dih,status='old')
3080 c do irec=1,maxres-3 ! loop for reading sigma_dih
3081 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3082 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3083 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3084 c & sigma_dih(k,i+nnt-1)
3089 if (idomain(k,i).eq.0) then
3093 dih(k,i)=phiref(i) ! right?
3094 c read (ientin,*) sigma_dih(k,i) ! original variant
3095 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3096 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3097 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3098 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3099 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3101 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3102 & rescore(k,i-2)+rescore(k,i-3))/4.0
3103 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3104 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3105 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3106 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3107 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3108 c Instead of res sim other local measure of b/b str reliability possible
3109 if (sigma_dih(k,i).ne.0)
3110 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3111 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3116 if (waga_theta.gt.0.0d0) then
3117 c open (ientin,file=tpl_k_sigma_theta,status='old')
3118 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3119 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3120 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3121 c & sigma_theta(k,i+nnt-1)
3126 do i = nnt+2,nct ! right? without parallel.
3127 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3128 c do i=ithet_start,ithet_end ! with FG parallel.
3129 if (idomain(k,i).eq.0) then
3130 sigma_theta(k,i)=0.0
3133 thetatpl(k,i)=thetaref(i)
3134 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3135 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3136 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3137 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3138 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3139 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3140 & rescore(k,i-2))/3.0
3141 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3142 if (sigma_theta(k,i).ne.0)
3143 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3145 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3146 c rescore(k,i-2) ! right expression ?
3147 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3151 if (waga_d.gt.0.0d0) then
3152 c open (ientin,file=tpl_k_sigma_d,status='old')
3153 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3154 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3155 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3156 c & sigma_d(k,i+nnt-1)
3160 do i = nnt,nct ! right? without parallel.
3161 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3162 c do i=loc_start,loc_end ! with FG parallel.
3163 if (itype(i).eq.10) cycle
3164 if (idomain(k,i).eq.0 ) then
3171 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3172 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3173 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3174 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3175 sigma_d(k,i)=rescore(k,i) ! right expression ?
3176 if (sigma_d(k,i).ne.0)
3177 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3179 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3180 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3181 c read (ientin,*) sigma_d(k,i) ! 1st variant
3182 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
3187 c remove distance restraints not used in any model from the list
3188 c shift data in all arrays
3190 if (waga_dist.ne.0.0d0) then
3196 if (ii_in_use(ii).eq.0.and.liiflag) then
3200 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3201 & .not.liiflag.and.ii.eq.lim_odl) then
3202 if (ii.eq.lim_odl) then
3203 iishift=ii-iistart+1
3208 do ki=iistart,lim_odl-iishift
3209 ires_homo(ki)=ires_homo(ki+iishift)
3210 jres_homo(ki)=jres_homo(ki+iishift)
3211 ii_in_use(ki)=ii_in_use(ki+iishift)
3212 do k=1,constr_homology
3213 odl(k,ki)=odl(k,ki+iishift)
3214 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3215 l_homo(k,ki)=l_homo(k,ki+iishift)
3219 lim_odl=lim_odl-iishift
3225 endif ! .not. klapaucjusz
3227 if (constr_homology.gt.0) call homology_partition
3228 if (constr_homology.gt.0) call init_int_table
3229 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
3230 cd & "lim_xx=",lim_xx
3231 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3232 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3236 if (.not.lprn) return
3237 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3238 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3239 write (iout,*) "Distance restraints from templates"
3241 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3242 & ii,ires_homo(ii),jres_homo(ii),
3243 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3244 & ki=1,constr_homology)
3246 write (iout,*) "Dihedral angle restraints from templates"
3248 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3249 & (rad2deg*dih(ki,i),
3250 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3252 write (iout,*) "Virtual-bond angle restraints from templates"
3254 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3255 & (rad2deg*thetatpl(ki,i),
3256 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3258 write (iout,*) "SC restraints from templates"
3260 write(iout,'(i5,100(4f8.2,4x))') i,
3261 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3262 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3265 c -----------------------------------------------------------------
3268 c -----------------------------------------------------------------
3269 subroutine read_klapaucjusz
3271 include 'DIMENSIONS'
3275 include 'COMMON.SETUP'
3276 include 'COMMON.CONTROL'
3277 include 'COMMON.CHAIN'
3278 include 'COMMON.IOUNITS'
3280 include 'COMMON.GEO'
3281 include 'COMMON.INTERACT'
3282 include 'COMMON.NAMES'
3283 character*256 fragfile
3284 integer ninclust(maxclust),inclust(max_template,maxclust),
3285 & nresclust(maxclust),iresclust(maxres,maxclust)
3288 character*24 model_ki_dist, model_ki_angle
3289 character*500 controlcard
3290 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3291 logical lprn /.true./
3297 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3298 double precision, dimension (max_template,maxres) :: rescore
3299 double precision, dimension (max_template,maxres) :: rescore2
3300 character*24 pdbfile,tpl_k_rescore
3303 c For new homol impl
3305 include 'COMMON.VAR'
3307 call getenv("FRAGFILE",fragfile)
3308 open(ientin,file=fragfile,status="old",err=10)
3309 read(ientin,*) constr_homology,nclust
3315 do k=1,constr_homology
3316 read(ientin,'(a)') pdbfile
3317 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3318 & pdbfile(:ilen(pdbfile))
3319 open(ipdbin,file=pdbfile,status='old',err=33)
3321 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3322 & pdbfile(:ilen(pdbfile))
3326 call readpdb_template(k)
3334 read(ientin,*) ninclust(i),nresclust(i)
3335 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3336 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3339 c Loop over clusters
3342 do ll = 1,ninclust(l)
3350 idomain(k,iresclust(i,l)+1) = 1
3352 idomain(k,iresclust(i,l)) = 1
3356 c Distance restraints
3359 C Copy the coordinates from reference coordinates (?)
3363 c write (iout,*) "c(",j,i,") =",c(j,i)
3366 call int_from_cart(.true.,.false.)
3367 call sc_loc_geom(.false.)
3369 thetaref(i)=theta(i)
3372 if (waga_dist.ne.0.0d0) then
3380 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3381 c write (iout,*) k,i,j,distal,dist2_cut
3383 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3384 & .and. distal.le.dist2_cut ) then
3390 c write (iout,*) "k",k
3391 c write (iout,*) "i",i," j",j," constr_homology",
3396 if (read2sigma) then
3399 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3401 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3402 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3403 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3405 if (odl(k,ii).le.dist_cut) then
3406 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3409 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3410 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3412 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3413 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3417 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3420 c l_homo(k,ii)=.false.
3427 c Theta, dihedral and SC retraints
3429 if (waga_angle.gt.0.0d0) then
3431 if (idomain(k,i).eq.0) then
3432 c sigma_dih(k,i)=0.0
3436 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3437 & rescore(k,i-2)+rescore(k,i-3))/4.0
3438 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3439 c & " sigma_dihed",sigma_dih(k,i)
3440 if (sigma_dih(k,i).ne.0)
3441 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3446 if (waga_theta.gt.0.0d0) then
3448 if (idomain(k,i).eq.0) then
3449 c sigma_theta(k,i)=0.0
3452 thetatpl(k,i)=thetaref(i)
3453 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3454 & rescore(k,i-2))/3.0
3455 if (sigma_theta(k,i).ne.0)
3456 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3460 if (waga_d.gt.0.0d0) then
3462 if (itype(i).eq.10) cycle
3463 if (idomain(k,i).eq.0 ) then
3470 sigma_d(k,i)=rescore(k,i)
3471 if (sigma_d(k,i).ne.0)
3472 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3473 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3479 c remove distance restraints not used in any model from the list
3480 c shift data in all arrays
3482 if (waga_dist.ne.0.0d0) then
3488 if (ii_in_use(ii).eq.0.and.liiflag) then
3492 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3493 & .not.liiflag.and.ii.eq.lim_odl) then
3494 if (ii.eq.lim_odl) then
3495 iishift=ii-iistart+1
3500 do ki=iistart,lim_odl-iishift
3501 ires_homo(ki)=ires_homo(ki+iishift)
3502 jres_homo(ki)=jres_homo(ki+iishift)
3503 ii_in_use(ki)=ii_in_use(ki+iishift)
3504 do k=1,constr_homology
3505 odl(k,ki)=odl(k,ki+iishift)
3506 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3507 l_homo(k,ki)=l_homo(k,ki+iishift)
3511 lim_odl=lim_odl-iishift
3518 10 stop "Error infragment file"
3520 c----------------------------------------------------------------------
3523 subroutine flush(iu)
3528 subroutine flush(iu)
3534 c------------------------------------------------------------------------------
3535 subroutine copy_to_tmp(source)
3536 include "DIMENSIONS"
3537 include "COMMON.IOUNITS"
3538 character*(*) source
3539 character* 256 tmpfile
3543 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3544 inquire(file=tmpfile,exist=ex)
3546 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3547 & " to temporary directory..."
3548 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3549 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3553 c------------------------------------------------------------------------------
3554 subroutine move_from_tmp(source)
3555 include "DIMENSIONS"
3556 include "COMMON.IOUNITS"
3557 character*(*) source
3560 write (*,*) "Moving ",source(:ilen(source)),
3561 & " from temporary directory to working directory"
3562 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3563 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3566 c------------------------------------------------------------------------------
3567 subroutine random_init(seed)
3569 C Initialize random number generator
3571 implicit real*8 (a-h,o-z)
3572 include 'DIMENSIONS'
3578 logical OKRandom, prng_restart
3580 integer iseed_array(4)
3582 include 'COMMON.IOUNITS'
3583 include 'COMMON.TIME1'
3584 include 'COMMON.THREAD'
3585 include 'COMMON.SBRIDGE'
3586 include 'COMMON.CONTROL'
3587 include 'COMMON.MCM'
3588 include 'COMMON.MAP'
3589 include 'COMMON.HEADER'
3590 csa include 'COMMON.CSA'
3591 include 'COMMON.CHAIN'
3592 include 'COMMON.MUCA'
3594 include 'COMMON.FFIELD'
3595 include 'COMMON.SETUP'
3596 iseed=-dint(dabs(seed))
3597 if (iseed.eq.0) then
3598 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3599 & 'Random seed undefined. The program will stop.'
3600 write (*,'(/80(1h*)/20x,a/80(1h*))')
3601 & 'Random seed undefined. The program will stop.'
3603 call mpi_finalize(mpi_comm_world,ierr)
3605 stop 'Bad random seed.'
3608 if (fg_rank.eq.0) then
3612 if(me.eq.king .or. .not. out1file)
3613 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3614 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3615 OKRandom = prng_restart(me,iseedi8)
3618 tmp=65536.0d0**(4-i)
3619 iseed_array(i) = dint(seed/tmp)
3620 seed=seed-iseed_array(i)*tmp
3622 if(me.eq.king .or. .not. out1file)
3623 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3624 & (iseed_array(i),i=1,4)
3625 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3626 & (iseed_array(i),i=1,4)
3627 OKRandom = prng_restart(me,iseed_array)
3630 r1=ran_number(0.0D0,1.0D0)
3631 if(me.eq.king .or. .not. out1file)
3632 & write (iout,*) 'ran_num',r1
3633 if (r1.lt.0.0d0) OKRandom=.false.
3635 if (.not.OKRandom) then
3636 write (iout,*) 'PRNG IS NOT WORKING!!!'
3637 print *,'PRNG IS NOT WORKING!!!'
3640 call mpi_abort(mpi_comm_world,error_msg,ierr)
3643 write (iout,*) 'too many processors for parallel prng'
3644 write (*,*) 'too many processors for parallel prng'
3652 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)