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)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
49 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
54 c print *,"Processor",myrank," leaves READRTNS"
57 C-------------------------------------------------------------------------------
58 subroutine read_control
62 implicit real*8 (a-h,o-z)
66 logical OKRandom, prng_restart
69 include 'COMMON.IOUNITS'
70 include 'COMMON.TIME1'
71 include 'COMMON.THREAD'
72 include 'COMMON.SBRIDGE'
73 include 'COMMON.CONTROL'
76 include 'COMMON.HEADER'
77 csa include 'COMMON.CSA'
78 include 'COMMON.CHAIN'
81 include 'COMMON.FFIELD'
82 include 'COMMON.SETUP'
83 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
84 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
86 character*320 controlcard
91 read (INP,'(a)') titel
92 call card_concat(controlcard)
93 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
94 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
95 call reada(controlcard,'SEED',seed,0.0D0)
96 call random_init(seed)
97 C Set up the time limit (caution! The time must be input in minutes!)
98 read_cart=index(controlcard,'READ_CART').gt.0
99 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
100 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
101 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
102 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
103 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
104 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
105 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
106 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
107 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
108 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
109 call reada(controlcard,'DRMS',drms,0.1D0)
110 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
111 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
112 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
113 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
114 write (iout,'(a,f10.1)')'DRMS = ',drms
115 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
116 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
118 call readi(controlcard,'NZ_START',nz_start,0)
119 call readi(controlcard,'NZ_END',nz_end,0)
120 call readi(controlcard,'IZ_SC',iz_sc,0)
122 safety = 60.0d0*safety
125 call reada(controlcard,"T_BATH",t_bath,300.0d0)
126 minim=(index(controlcard,'MINIMIZE').gt.0)
127 dccart=(index(controlcard,'CART').gt.0)
128 overlapsc=(index(controlcard,'OVERLAP').gt.0)
129 overlapsc=.not.overlapsc
130 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
131 searchsc=.not.searchsc
132 sideadd=(index(controlcard,'SIDEADD').gt.0)
133 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
134 outpdb=(index(controlcard,'PDBOUT').gt.0)
135 outx=(index(controlcard,'XOUT').gt.0)
136 outmol2=(index(controlcard,'MOL2OUT').gt.0)
137 pdbref=(index(controlcard,'PDBREF').gt.0)
138 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
139 indpdb=index(controlcard,'PDBSTART')
140 extconf=(index(controlcard,'EXTCONF').gt.0)
141 call readi(controlcard,'IPRINT',iprint,0)
142 call readi(controlcard,'MAXGEN',maxgen,10000)
143 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
144 call readi(controlcard,"KDIAG",kdiag,0)
145 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
146 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
147 & write (iout,*) "RESCALE_MODE",rescale_mode
148 split_ene=index(controlcard,'SPLIT_ENE').gt.0
149 if (index(controlcard,'REGULAR').gt.0.0D0) then
150 call reada(controlcard,'WEIDIS',weidis,0.1D0)
154 call reada(controlcard,"CHECKGRAD_INC",checkgrad_inc,1.0d-4)
155 if (index(controlcard,'CHECKGRAD').gt.0) then
157 if (index(controlcard,'CART').gt.0) then
159 elseif (index(controlcard,'CARINT').gt.0) then
164 elseif (index(controlcard,'THREAD').gt.0) then
166 call readi(controlcard,'THREAD',nthread,0)
167 if (nthread.gt.0) then
168 call reada(controlcard,'WEIDIS',weidis,0.1D0)
171 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
172 stop 'Error termination in Read_Control.'
174 else if (index(controlcard,'MCMA').gt.0) then
176 else if (index(controlcard,'MCEE').gt.0) then
178 else if (index(controlcard,'MULTCONF').gt.0) then
180 else if (index(controlcard,'MAP').gt.0) then
182 call readi(controlcard,'MAP',nmap,0)
183 else if (index(controlcard,'CSA').gt.0) then
184 write(*,*) "CSA not supported in this version"
187 crc else if (index(controlcard,'ZSCORE').gt.0) then
189 crc ZSCORE is rm from UNRES, modecalc=9 is available
192 cfcm else if (index(controlcard,'MCMF').gt.0) then
194 else if (index(controlcard,'SOFTREG').gt.0) then
196 else if (index(controlcard,'CHECK_BOND').gt.0) then
198 else if (index(controlcard,'TEST').gt.0) then
200 else if (index(controlcard,'MD').gt.0) then
202 else if (index(controlcard,'RE ').gt.0) then
206 lmuca=index(controlcard,'MUCA').gt.0
207 call readi(controlcard,'MUCADYN',mucadyn,0)
208 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
209 if (lmuca .and. (me.eq.king .or. .not.out1file ))
211 write (iout,*) 'MUCADYN=',mucadyn
212 write (iout,*) 'MUCASMOOTH=',muca_smooth
215 iscode=index(controlcard,'ONE_LETTER')
216 indphi=index(controlcard,'PHI')
217 indback=index(controlcard,'BACK')
218 iranconf=index(controlcard,'RAND_CONF')
219 i2ndstr=index(controlcard,'USE_SEC_PRED')
220 gradout=index(controlcard,'GRADOUT').gt.0
221 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
223 if(me.eq.king.or..not.out1file)
224 & write (iout,'(2a)') diagmeth(kdiag),
225 & ' routine used to diagonalize matrices.'
228 c--------------------------------------------------------------------------
229 subroutine read_REMDpar
233 implicit real*8 (a-h,o-z)
235 include 'COMMON.IOUNITS'
236 include 'COMMON.TIME1'
239 include 'COMMON.LANGEVIN'
241 include 'COMMON.LANGEVIN.lang0'
243 include 'COMMON.INTERACT'
244 include 'COMMON.NAMES'
246 include 'COMMON.REMD'
247 include 'COMMON.CONTROL'
248 include 'COMMON.SETUP'
250 character*320 controlcard
251 character*3200 controlcard1
252 integer iremd_m_total
254 if(me.eq.king.or..not.out1file)
255 & write (iout,*) "REMD setup"
257 call card_concat(controlcard)
258 call readi(controlcard,"NREP",nrep,3)
259 call readi(controlcard,"NSTEX",nstex,1000)
260 call reada(controlcard,"RETMIN",retmin,10.0d0)
261 call reada(controlcard,"RETMAX",retmax,1000.0d0)
262 mremdsync=(index(controlcard,'SYNC').gt.0)
263 call readi(controlcard,"NSYN",i_sync_step,100)
264 restart1file=(index(controlcard,'REST1FILE').gt.0)
265 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
266 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
267 if(max_cache_traj_use.gt.max_cache_traj)
268 & max_cache_traj_use=max_cache_traj
269 if(me.eq.king.or..not.out1file) then
270 cd if (traj1file) then
271 crc caching is in testing - NTWX is not ignored
272 cd write (iout,*) "NTWX value is ignored"
273 cd write (iout,*) " trajectory is stored to one file by master"
274 cd write (iout,*) " before exchange at NSTEX intervals"
276 write (iout,*) "NREP= ",nrep
277 write (iout,*) "NSTEX= ",nstex
278 write (iout,*) "SYNC= ",mremdsync
279 write (iout,*) "NSYN= ",i_sync_step
280 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
283 t_exchange_only=(index(controlcard,'TONLY').gt.0)
284 call readi(controlcard,"HREMD",hremd,0)
285 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
286 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
288 if(usampl.and.hremd.gt.0) then
290 & "========== ERROR: USAMPL and HREMD cannot be used together"
292 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
299 if (index(controlcard,'TLIST').gt.0) then
301 call card_concat(controlcard1)
302 read(controlcard1,*) (remd_t(i),i=1,nrep)
303 if(me.eq.king.or..not.out1file)
304 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
307 if (index(controlcard,'MLIST').gt.0) then
309 call card_concat(controlcard1)
310 read(controlcard1,*) (remd_m(i),i=1,nrep)
311 if(me.eq.king.or..not.out1file) then
312 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
315 iremd_m_total=iremd_m_total+remd_m(i)
318 write (iout,*) 'Total number of replicas ',
319 & iremd_m_total*hremd
321 write (iout,*) 'Total number of replicas ',iremd_m_total
325 if(me.eq.king.or..not.out1file)
326 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
329 c--------------------------------------------------------------------------
330 subroutine read_MDpar
334 implicit real*8 (a-h,o-z)
336 include 'COMMON.IOUNITS'
337 include 'COMMON.TIME1'
340 include 'COMMON.LANGEVIN'
342 include 'COMMON.LANGEVIN.lang0'
344 include 'COMMON.INTERACT'
345 include 'COMMON.NAMES'
347 include 'COMMON.SETUP'
348 include 'COMMON.CONTROL'
349 include 'COMMON.SPLITELE'
351 character*320 controlcard
353 call card_concat(controlcard)
354 call readi(controlcard,"NSTEP",n_timestep,1000000)
355 call readi(controlcard,"NTWE",ntwe,100)
356 call readi(controlcard,"NTWX",ntwx,1000)
357 call reada(controlcard,"DT",d_time,1.0d-1)
358 call reada(controlcard,"DVMAX",dvmax,2.0d1)
359 call reada(controlcard,"DAMAX",damax,1.0d1)
360 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
361 call readi(controlcard,"LANG",lang,0)
362 RESPA = index(controlcard,"RESPA") .gt. 0
363 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
364 ntime_split0=ntime_split
365 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
366 ntime_split0=ntime_split
367 call reada(controlcard,"R_CUT",r_cut,2.0d0)
368 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
369 rest = index(controlcard,"REST").gt.0
370 tbf = index(controlcard,"TBF").gt.0
371 call readi(controlcard,"HMC",hmc,0)
372 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
373 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
374 tnh = index(controlcard,"NOSEHOOVER96").gt.0
375 if (RESPA.and.tnh)then
376 xiresp = index(controlcard,"XIRESP").gt.0
378 call reada(controlcard,"Q_NP",Q_np,0.1d0)
379 usampl = index(controlcard,"USAMPL").gt.0
380 mdpdb = index(controlcard,"MDPDB").gt.0
381 call reada(controlcard,"T_BATH",t_bath,300.0d0)
382 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
383 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
384 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
385 if (count_reset_moment.eq.0) count_reset_moment=1000000000
386 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
387 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
388 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
389 if (count_reset_vel.eq.0) count_reset_vel=1000000000
390 large = index(controlcard,"LARGE").gt.0
391 print_compon = index(controlcard,"PRINT_COMPON").gt.0
392 rattle = index(controlcard,"RATTLE").gt.0
393 preminim = index(controlcard,"PREMINIM").gt.0
395 dccart=(index(controlcard,'CART').gt.0)
398 c if performing umbrella sampling, fragments constrained are read from the fragment file
404 if(me.eq.king.or..not.out1file) then
406 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
408 write (iout,'(a)') "The units are:"
409 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
410 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
411 & " acceleration: angstrom/(48.9 fs)**2"
412 write (iout,'(a)') "energy: kcal/mol, temperature: K"
414 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
415 write (iout,'(a60,f10.5,a)')
416 & "Initial time step of numerical integration:",d_time,
418 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
420 write (iout,'(2a,i4,a)')
421 & "A-MTS algorithm used; initial time step for fast-varying",
422 & " short-range forces split into",ntime_split," steps."
423 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
424 & r_cut," lambda",rlamb
426 write (iout,'(2a,f10.5)')
427 & "Maximum acceleration threshold to reduce the time step",
428 & "/increase split number:",damax
429 write (iout,'(2a,f10.5)')
430 & "Maximum predicted energy drift to reduce the timestep",
431 & "/increase split number:",edriftmax
432 write (iout,'(a60,f10.5)')
433 & "Maximum velocity threshold to reduce velocities:",dvmax
434 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
435 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
436 if (rattle) write (iout,'(a60)')
437 & "Rattle algorithm used to constrain the virtual bonds"
438 if (preminim .or. iranconf.gt.0) then
440 & "Initial structure will be energy-minimized"
445 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
446 call reada(controlcard,"RWAT",rwat,1.4d0)
447 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
448 surfarea=index(controlcard,"SURFAREA").gt.0
449 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
450 if(me.eq.king.or..not.out1file)then
451 write (iout,'(/a,$)') "Langevin dynamics calculation"
454 & " with direct integration of Langevin equations"
455 else if (lang.eq.2) then
456 write (iout,'(a/)') " with TINKER stochasic MD integrator"
457 else if (lang.eq.3) then
458 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
459 else if (lang.eq.4) then
460 write (iout,'(a/)') " in overdamped mode"
462 write (iout,'(//a,i5)')
463 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
466 write (iout,'(a60,f10.5)') "Temperature:",t_bath
467 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
468 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
469 write (iout,'(a60,f10.5)')
470 & "Scaling factor of the friction forces:",scal_fric
471 if (surfarea) write (iout,'(2a,i10,a)')
472 & "Friction coefficients will be scaled by solvent-accessible",
473 & " surface area every",reset_fricmat," steps."
475 c Calculate friction coefficients and bounds of stochastic forces
476 eta=6*pi*cPoise*etawat
477 if(me.eq.king.or..not.out1file)
478 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
480 gamp=scal_fric*(pstok+rwat)*eta
481 stdfp=dsqrt(2*Rb*t_bath/d_time)
483 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
484 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
486 if(me.eq.king.or..not.out1file)then
487 write (iout,'(/2a/)')
488 & "Radii of site types and friction coefficients and std's of",
489 & " stochastic forces of fully exposed sites"
490 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
492 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
493 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
497 if(me.eq.king.or..not.out1file)then
498 write (iout,'(a)') "Berendsen bath calculation"
499 write (iout,'(a60,f10.5)') "Temperature:",t_bath
500 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
502 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
503 & count_reset_moment," steps"
505 & write (iout,'(a,i10,a)')
506 & "Velocities will be reset at random every",count_reset_vel,
509 else if (tnp .or. tnp1 .or. tnh) then
510 if (tnp .or. tnp1) then
511 write (iout,'(a)') "Nose-Poincare bath calculation"
512 if (tnp) write (iout,'(a)')
513 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
514 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
516 write (iout,'(a)') "Nose-Hoover bath calculation"
517 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
527 WDTI(i) = 1.0*d_time/nresn
534 write (iout,'(a)') "NVT-XI-RESPA algorithm"
536 write (iout,'(a)') "NVT-XO-RESPA algorithm"
539 WDTIi(i) = 1.0*d_time/nresn/ntime_split
547 write (iout,'(a60,f10.5)') "Temperature:",t_bath
548 write (iout,'(a60,f10.5)') "Q =",Q_np
550 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
551 & count_reset_moment," steps"
553 & write (iout,'(a,i10,a)')
554 & "Velocities will be reset at random every",count_reset_vel,
557 else if (hmc.gt.0) then
558 write (iout,'(a)') "Hybrid Monte Carlo calculation"
559 write (iout,'(a60,f10.5)') "Temperature:",t_bath
560 write (iout,'(a60,i10)')
561 & "Number of MD steps between Metropolis tests:",hmc
564 if(me.eq.king.or..not.out1file)
565 & write (iout,'(a31)') "Microcanonical mode calculation"
567 if(me.eq.king.or..not.out1file)then
568 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
570 write(iout,*) "MD running with constraints."
571 write(iout,*) "Equilibration time ", eq_time, " mtus."
572 write(iout,*) "Constraining ", nfrag," fragments."
573 write(iout,*) "Length of each fragment, weight and q0:"
575 write (iout,*) "Set of restraints #",iset
577 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
578 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
580 write(iout,*) "constraints between ", npair, "fragments."
581 write(iout,*) "constraint pairs, weights and q0:"
583 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
584 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
586 write(iout,*) "angle constraints within ", nfrag_back,
587 & "backbone fragments."
588 write(iout,*) "fragment, weights:"
590 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
591 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
592 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
595 iset=mod(kolor,nset)+1
598 if(me.eq.king.or..not.out1file)
599 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
602 c------------------------------------------------------------------------------
605 C Read molecular data.
607 implicit real*8 (a-h,o-z)
613 include 'COMMON.IOUNITS'
616 include 'COMMON.INTERACT'
617 include 'COMMON.LOCAL'
618 include 'COMMON.NAMES'
619 include 'COMMON.CHAIN'
620 include 'COMMON.FFIELD'
621 include 'COMMON.SBRIDGE'
622 include 'COMMON.HEADER'
623 include 'COMMON.CONTROL'
624 include 'COMMON.DBASE'
625 include 'COMMON.THREAD'
626 include 'COMMON.CONTACTS'
627 include 'COMMON.TORCNSTR'
628 include 'COMMON.TIME1'
629 include 'COMMON.BOUNDS'
631 include 'COMMON.REMD'
632 include 'COMMON.SETUP'
633 character*4 sequence(maxres)
635 double precision x(maxvar)
636 character*256 pdbfile
637 character*320 weightcard
638 character*80 weightcard_t,ucase
639 dimension itype_pdb(maxres)
640 common /pizda/ itype_pdb
641 logical seq_comp,fail
642 double precision energia(0:n_ene)
649 C Read weights of the subsequent energy terms.
662 if(me.eq.king.or..not.out1file) then
663 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
664 write (iout,*) 'Current weights for processor ',
665 & me,' set ',i2set(me)
669 call card_concat(weightcard)
670 call reada(weightcard,'WLONG',wlong,1.0D0)
671 call reada(weightcard,'WSC',wsc,wlong)
672 call reada(weightcard,'WSCP',wscp,wlong)
673 call reada(weightcard,'WELEC',welec,1.0D0)
674 call reada(weightcard,'WVDWPP',wvdwpp,welec)
675 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
676 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
677 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
678 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
679 call reada(weightcard,'WTURN3',wturn3,1.0D0)
680 call reada(weightcard,'WTURN4',wturn4,1.0D0)
681 call reada(weightcard,'WTURN6',wturn6,1.0D0)
682 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
683 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
684 call reada(weightcard,'WBOND',wbond,1.0D0)
685 call reada(weightcard,'WTOR',wtor,1.0D0)
686 call reada(weightcard,'WTORD',wtor_d,1.0D0)
687 call reada(weightcard,'WANG',wang,1.0D0)
688 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
689 call reada(weightcard,'SCAL14',scal14,0.4D0)
690 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
691 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
692 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
693 call reada(weightcard,'TEMP0',temp0,300.0d0)
694 if (index(weightcard,'SOFT').gt.0) ipot=6
695 C 12/1/95 Added weight for the multi-body term WCORR
696 call reada(weightcard,'WCORRH',wcorr,1.0D0)
697 if (wcorr4.gt.0.0d0) wcorr=wcorr4
705 hweights(i,7)=wel_loc
708 hweights(i,10)=wturn6
710 hweights(i,12)=wscloc
712 hweights(i,14)=wtor_d
713 hweights(i,15)=wstrain
714 hweights(i,16)=wvdwpp
716 hweights(i,18)=scal14
717 hweights(i,21)=wsccor
722 weights(i)=hweights(i2set(me),i)
746 call card_concat(weightcard)
747 call reada(weightcard,'WLONG',wlong,1.0D0)
748 call reada(weightcard,'WSC',wsc,wlong)
749 call reada(weightcard,'WSCP',wscp,wlong)
750 call reada(weightcard,'WELEC',welec,1.0D0)
751 call reada(weightcard,'WVDWPP',wvdwpp,welec)
752 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
753 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
754 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
755 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
756 call reada(weightcard,'WTURN3',wturn3,1.0D0)
757 call reada(weightcard,'WTURN4',wturn4,1.0D0)
758 call reada(weightcard,'WTURN6',wturn6,1.0D0)
759 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
760 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
761 call reada(weightcard,'WBOND',wbond,1.0D0)
762 call reada(weightcard,'WTOR',wtor,1.0D0)
763 call reada(weightcard,'WTORD',wtor_d,1.0D0)
764 call reada(weightcard,'WANG',wang,1.0D0)
765 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
766 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
767 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
768 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
769 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
770 call reada(weightcard,'SCAL14',scal14,0.4D0)
771 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
772 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
773 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
774 call reada(weightcard,'TEMP0',temp0,300.0d0)
775 if (index(weightcard,'SOFT').gt.0) ipot=6
776 C 12/1/95 Added weight for the multi-body term WCORR
777 call reada(weightcard,'WCORRH',wcorr,1.0D0)
778 if (wcorr4.gt.0.0d0) wcorr=wcorr4
799 weights(25)=wdfa_dist
802 weights(28)=wdfa_beta
804 if(me.eq.king.or..not.out1file)
805 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
806 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
808 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
809 10 format (/'Energy-term weights (unscaled):'//
810 & 'WSCC= ',f10.6,' (SC-SC)'/
811 & 'WSCP= ',f10.6,' (SC-p)'/
812 & 'WELEC= ',f10.6,' (p-p electr)'/
813 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
814 & 'WBOND= ',f10.6,' (stretching)'/
815 & 'WANG= ',f10.6,' (bending)'/
816 & 'WSCLOC= ',f10.6,' (SC local)'/
817 & 'WTOR= ',f10.6,' (torsional)'/
818 & 'WTORD= ',f10.6,' (double torsional)'/
819 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
820 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
821 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
822 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
823 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
824 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
825 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
826 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
827 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
828 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
829 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
830 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
831 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
832 if(me.eq.king.or..not.out1file)then
833 if (wcorr4.gt.0.0d0) then
834 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
835 & 'between contact pairs of peptide groups'
836 write (iout,'(2(a,f5.3/))')
837 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
838 & 'Range of quenching the correlation terms:',2*delt_corr
839 else if (wcorr.gt.0.0d0) then
840 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
841 & 'between contact pairs of peptide groups'
843 write (iout,'(a,f8.3)')
844 & 'Scaling factor of 1,4 SC-p interactions:',scal14
845 write (iout,'(a,f8.3)')
846 & 'General scaling factor of SC-p interactions:',scalscp
848 r0_corr=cutoff_corr-delt_corr
850 aad(i,1)=scalscp*aad(i,1)
851 aad(i,2)=scalscp*aad(i,2)
852 bad(i,1)=scalscp*bad(i,1)
853 bad(i,2)=scalscp*bad(i,2)
855 call rescale_weights(t_bath)
856 if(me.eq.king.or..not.out1file)
857 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
858 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
860 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
861 22 format (/'Energy-term weights (scaled):'//
862 & 'WSCC= ',f10.6,' (SC-SC)'/
863 & 'WSCP= ',f10.6,' (SC-p)'/
864 & 'WELEC= ',f10.6,' (p-p electr)'/
865 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
866 & 'WBOND= ',f10.6,' (stretching)'/
867 & 'WANG= ',f10.6,' (bending)'/
868 & 'WSCLOC= ',f10.6,' (SC local)'/
869 & 'WTOR= ',f10.6,' (torsional)'/
870 & 'WTORD= ',f10.6,' (double torsional)'/
871 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
872 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
873 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
874 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
875 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
876 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
877 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
878 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
879 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
880 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
881 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
882 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
883 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
884 if(me.eq.king.or..not.out1file)
885 & write (iout,*) "Reference temperature for weights calculation:",
887 call reada(weightcard,"D0CM",d0cm,3.78d0)
888 call reada(weightcard,"AKCM",akcm,15.1d0)
889 call reada(weightcard,"AKTH",akth,11.0d0)
890 call reada(weightcard,"AKCT",akct,12.0d0)
891 call reada(weightcard,"V1SS",v1ss,-1.08d0)
892 call reada(weightcard,"V2SS",v2ss,7.61d0)
893 call reada(weightcard,"V3SS",v3ss,13.7d0)
894 call reada(weightcard,"EBR",ebr,-5.50D0)
895 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
897 dyn_ss_mask(i)=.false.
901 dyn_ssbond_ij(i,j)=1.0d300
904 call reada(weightcard,"HT",Ht,0.0D0)
906 ss_depth=ebr/wsc-0.25*eps(1,1)
907 Ht=Ht/wsc-0.25*eps(1,1)
908 akcm=akcm*wstrain/wsc
909 akth=akth*wstrain/wsc
910 akct=akct*wstrain/wsc
911 v1ss=v1ss*wstrain/wsc
912 v2ss=v2ss*wstrain/wsc
913 v3ss=v3ss*wstrain/wsc
915 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
918 if(me.eq.king.or..not.out1file) then
919 write (iout,*) "Parameters of the SS-bond potential:"
920 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
922 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
923 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
924 write (iout,*)" HT",Ht
925 print *,'indpdb=',indpdb,' pdbref=',pdbref
927 if (indpdb.gt.0 .or. pdbref) then
928 read(inp,'(a)') pdbfile
929 if(me.eq.king.or..not.out1file)
930 & write (iout,'(2a)') 'PDB data will be read from file ',
931 & pdbfile(:ilen(pdbfile))
932 open(ipdbin,file=pdbfile,status='old',err=33)
934 33 write (iout,'(a)') 'Error opening PDB file.'
937 c print *,'Begin reading pdb data'
946 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
947 & (crefjlee(j,i+nres),j=1,3)
950 c print *,'Finished reading pdb data'
951 if(me.eq.king.or..not.out1file)
952 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
953 & ' nstart_sup=',nstart_sup
955 itype_pdb(i)=itype(i)
959 nct=nstart_sup+nsup-1
960 call contact(.false.,ncont_ref,icont_ref,co)
963 C Following 2 lines for diagnostics; comment out if not needed
964 write (iout,*) "Before sideadd"
966 if(me.eq.king.or..not.out1file)
967 & write(iout,*)'Adding sidechains'
974 do while (fail.and.nsi.le.maxsi)
975 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
978 if(fail) write(iout,*)'Adding sidechain failed for res ',
979 & i,' after ',nsi,' trials'
982 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
985 C Following 2 lines for diagnostics; comment out if not needed
986 c write (iout,*) "After sideadd"
989 if (indpdb.eq.0) then
990 C Read sequence if not taken from the pdb file.
992 c print *,'nres=',nres
993 if (iscode.gt.0) then
994 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
996 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
998 C Convert sequence to numeric code
1000 itype(i)=rescode(i,sequence(i),iscode)
1002 if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
1004 & "Glycine is the first full residue, initial dummy deleted"
1010 if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
1012 & "Glycine is the last full residue, terminal dummy deleted"
1016 C Assign initial virtual bond lengths
1022 vbld(i+nres)=dsc(itype(i))
1023 vbld_inv(i+nres)=dsc_inv(itype(i))
1024 c write (iout,*) "i",i," itype",itype(i),
1025 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
1029 c print '(20i4)',(itype(i),i=1,nres)
1032 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
1034 if (itype(i).eq.21) then
1038 else if (itype(i+1).ne.20) then
1040 else if (itype(i).ne.20) then
1047 if(me.eq.king.or..not.out1file)then
1048 write (iout,*) "ITEL"
1050 write (iout,*) i,itype(i),itel(i)
1052 print *,'Call Read_Bridge.'
1055 C 8/13/98 Set limits to generating the dihedral angles
1060 read (inp,*) ndih_constr
1061 if (ndih_constr.gt.0) then
1063 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1064 if(me.eq.king.or..not.out1file)then
1066 & 'There are',ndih_constr,' constraints on phi angles.'
1068 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1072 phi0(i)=deg2rad*phi0(i)
1073 drange(i)=deg2rad*drange(i)
1075 if(me.eq.king.or..not.out1file)
1076 & write (iout,*) 'FTORS',ftors
1079 phibound(1,ii) = phi0(i)-drange(i)
1080 phibound(2,ii) = phi0(i)+drange(i)
1085 if (me.eq.king) then
1087 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1089 write (iout,'(a3,i5,2f10.1)')
1090 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1096 cd print *,'NNT=',NNT,' NCT=',NCT
1097 if (itype(1).eq.21) nnt=2
1098 if (itype(nres).eq.21) nct=nct-1
1100 C Bartek:READ init_vars
1101 C Initialize variables!
1102 C Juyong:READ read_info
1103 C READ fragment information!!
1104 C both routines should be in dfa.F file!!
1107 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1108 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1110 print*, 'init_dfa_vars finished!'
1112 print*, 'read_dfa_info finished!'
1117 if(me.eq.king.or..not.out1file)
1118 & write (iout,'(a,i3)') 'nsup=',nsup
1120 if (nsup.le.(nct-nnt+1)) then
1121 do i=0,nct-nnt+1-nsup
1122 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1128 & 'Error - sequences to be superposed do not match.'
1131 do i=0,nsup-(nct-nnt+1)
1132 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1134 nstart_sup=nstart_sup+i
1140 & 'Error - sequences to be superposed do not match.'
1143 if (nsup.eq.0) nsup=nct-nnt
1144 if (nstart_sup.eq.0) nstart_sup=nnt
1145 if (nstart_seq.eq.0) nstart_seq=nnt
1146 if(me.eq.king.or..not.out1file)
1147 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1148 & ' nstart_seq=',nstart_seq
1150 c--- Zscore rms -------
1151 if (nz_start.eq.0) nz_start=nnt
1152 if (nz_end.eq.0 .and. nsup.gt.0) then
1154 else if (nz_end.eq.0) then
1157 if(me.eq.king.or..not.out1file)then
1158 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1159 write (iout,*) 'IZ_SC=',iz_sc
1161 c----------------------
1164 if (.not.pdbref) then
1165 call read_angles(inp,*38)
1167 38 write (iout,'(a)') 'Error reading reference structure.'
1169 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1170 stop 'Error reading reference structure'
1174 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1183 call contact(.true.,ncont_ref,icont_ref,co)
1185 if(me.eq.king.or..not.out1file)
1186 & write (iout,*) 'Contact order:',co
1188 if(me.eq.king.or..not.out1file)
1189 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1192 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1194 if(me.eq.king.or..not.out1file)
1195 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1196 & icont_ref(1,i),' ',
1197 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1201 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1202 if (constr_dist.gt.0) then
1203 call read_dist_constr
1207 if (constr_homology.gt.0) then
1208 call read_constr_homology
1209 if (indpdb.gt.0 .or. pdbref) then
1212 c(j,i)=crefjlee(j,i)
1213 cref(j,i)=crefjlee(j,i)
1218 write (iout,*) "Array C"
1220 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1221 & (c(j,i+nres),j=1,3)
1223 write (iout,*) "Array Cref"
1225 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1226 & (cref(j,i+nres),j=1,3)
1229 call int_from_cart1(.false.)
1230 call sc_loc_geom(.false.)
1232 thetaref(i)=theta(i)
1237 dc(j,i)=c(j,i+1)-c(j,i)
1238 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1243 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1244 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1252 if (nhpb.gt.0) call hpb_partition
1253 c write (iout,*) "After read_dist_constr nhpb",nhpb
1255 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1256 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1257 & modecalc.ne.10) then
1258 C If input structure hasn't been supplied from the PDB file read or generate
1260 if (iranconf.eq.0 .and. .not. extconf) then
1261 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1262 & write (iout,'(a)') 'Initial geometry will be read in.'
1264 read (inp,*) time,potE,uconst,t_bath,
1265 & nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn)
1266 read(inp,'(8f10.5)',end=36,err=36)
1267 & ((c(l,k),l=1,3),k=1,nres),
1268 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1269 call int_from_cart1(.false.)
1272 dc(j,i)=c(j,i+1)-c(j,i)
1273 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1277 if (itype(i).ne.10) then
1279 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1280 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1286 call read_angles(inp,*36)
1289 36 write (iout,'(a)') 'Error reading angle file.'
1291 call mpi_finalize( MPI_COMM_WORLD,IERR )
1293 stop 'Error reading angle file.'
1295 else if (extconf) then
1296 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1297 & write (iout,'(a)') 'Extended chain initial geometry.'
1299 theta(i)=90d0*deg2rad
1302 phi(i)=180d0*deg2rad
1305 alph(i)=110d0*deg2rad
1308 omeg(i)=-120d0*deg2rad
1311 if(me.eq.king.or..not.out1file)
1312 & write (iout,'(a)') 'Random-generated initial geometry.'
1316 if (me.eq.king .or. fg_rank.eq.0 .and. (
1317 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1321 call gen_rand_conf(itmp,*30)
1323 30 write (iout,*) 'Failed to generate random conformation',
1324 & ', itrial=',itrial
1325 write (*,*) 'Processor:',me,
1326 & ' Failed to generate random conformation',
1336 write (iout,'(a,i3,a)') 'Processor:',me,
1337 & ' error in generating random conformation.'
1338 write (*,'(a,i3,a)') 'Processor:',me,
1339 & ' error in generating random conformation.'
1342 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1349 elseif (modecalc.eq.4) then
1350 read (inp,'(a)') intinname
1351 open (intin,file=intinname,status='old',err=333)
1352 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1353 & write (iout,'(a)') 'intinname',intinname
1354 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1356 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1358 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1360 stop 'Error opening angle file.'
1364 C Generate distance constraints, if the PDB structure is to be regularized.
1365 if (nthread.gt.0) then
1366 call read_threadbase
1368 write (iout,*) "READRTNS: Calling setup_var"
1371 if (me.eq.king .or. .not. out1file)
1373 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1374 write (iout,'(/a,i3,a)')
1375 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1376 write (iout,'(20i4)') (iss(i),i=1,ns)
1378 write(iout,*)"Running with dynamic disulfide-bond formation"
1380 write (iout,'(/a/)') 'Pre-formed links are:'
1386 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1387 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1393 if (ns.gt.0.and.dyn_ss) then
1397 forcon(i-nss)=forcon(i)
1404 dyn_ss_mask(iss(i))=.true.
1407 if (i2ndstr.gt.0) call secstrp2dihc
1408 c call geom_to_var(nvar,x)
1409 c call etotal(energia(0))
1410 c call enerprint(energia(0))
1411 c call briefout(0,etot)
1413 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1414 cd write (iout,'(a)') 'Variable list:'
1415 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1417 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1418 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1419 & 'Processor',myrank,': end reading molecular data.'
1423 c--------------------------------------------------------------------------
1424 logical function seq_comp(itypea,itypeb,length)
1426 integer length,itypea(length),itypeb(length)
1429 if (itypea(i).ne.itypeb(i)) then
1437 c-----------------------------------------------------------------------------
1438 subroutine read_bridge
1439 C Read information about disulfide bridges.
1440 implicit real*8 (a-h,o-z)
1441 include 'DIMENSIONS'
1445 include 'COMMON.IOUNITS'
1446 include 'COMMON.GEO'
1447 include 'COMMON.VAR'
1448 include 'COMMON.INTERACT'
1449 include 'COMMON.LOCAL'
1450 include 'COMMON.NAMES'
1451 include 'COMMON.CHAIN'
1452 include 'COMMON.FFIELD'
1453 include 'COMMON.SBRIDGE'
1454 include 'COMMON.HEADER'
1455 include 'COMMON.CONTROL'
1456 include 'COMMON.DBASE'
1457 include 'COMMON.THREAD'
1458 include 'COMMON.TIME1'
1459 include 'COMMON.SETUP'
1460 C Read bridging residues.
1461 read (inp,*) ns,(iss(i),i=1,ns)
1463 if(me.eq.king.or..not.out1file)
1464 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1465 C Check whether the specified bridging residues are cystines.
1467 if (itype(iss(i)).ne.1) then
1468 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1469 & 'Do you REALLY think that the residue ',
1470 & restyp(itype(iss(i))),i,
1471 & ' can form a disulfide bridge?!!!'
1472 write (*,'(2a,i3,a)')
1473 & 'Do you REALLY think that the residue ',
1474 & restyp(itype(iss(i))),i,
1475 & ' can form a disulfide bridge?!!!'
1477 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1482 C Read preformed bridges.
1484 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1486 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1489 C Check if the residues involved in bridges are in the specified list of
1490 C bridging residues.
1493 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1494 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1495 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1496 & ' contains residues present in other pairs.'
1497 write (*,'(a,i3,a)') 'Disulfide pair',i,
1498 & ' contains residues present in other pairs.'
1500 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1506 if (ihpb(i).eq.iss(j)) goto 10
1508 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1511 if (jhpb(i).eq.iss(j)) goto 20
1513 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1519 ihpb(i)=ihpb(i)+nres
1520 jhpb(i)=jhpb(i)+nres
1526 c----------------------------------------------------------------------------
1527 subroutine read_x(kanal,*)
1528 implicit real*8 (a-h,o-z)
1529 include 'DIMENSIONS'
1530 include 'COMMON.GEO'
1531 include 'COMMON.VAR'
1532 include 'COMMON.CHAIN'
1533 include 'COMMON.IOUNITS'
1534 include 'COMMON.CONTROL'
1535 include 'COMMON.LOCAL'
1536 include 'COMMON.INTERACT'
1537 c Read coordinates from input
1539 read(kanal,'(8f10.5)',end=10,err=10)
1540 & ((c(l,k),l=1,3),k=1,nres),
1541 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1544 c(j,2*nres)=c(j,nres)
1546 call int_from_cart1(.false.)
1549 dc(j,i)=c(j,i+1)-c(j,i)
1550 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1554 if (itype(i).ne.10) then
1556 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1557 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1565 c----------------------------------------------------------------------------
1566 subroutine read_threadbase
1567 implicit real*8 (a-h,o-z)
1568 include 'DIMENSIONS'
1569 include 'COMMON.IOUNITS'
1570 include 'COMMON.GEO'
1571 include 'COMMON.VAR'
1572 include 'COMMON.INTERACT'
1573 include 'COMMON.LOCAL'
1574 include 'COMMON.NAMES'
1575 include 'COMMON.CHAIN'
1576 include 'COMMON.FFIELD'
1577 include 'COMMON.SBRIDGE'
1578 include 'COMMON.HEADER'
1579 include 'COMMON.CONTROL'
1580 include 'COMMON.DBASE'
1581 include 'COMMON.THREAD'
1582 include 'COMMON.TIME1'
1583 C Read pattern database for threading.
1584 read (icbase,*) nseq
1586 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1587 & nres_base(2,i),nres_base(3,i)
1588 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1590 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1591 c & nres_base(2,i),nres_base(3,i)
1592 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1596 if (weidis.eq.0.0D0) weidis=0.1D0
1605 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1606 write (iout,'(a,i5)') 'nexcl: ',nexcl
1607 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1610 c------------------------------------------------------------------------------
1611 subroutine setup_var
1612 implicit real*8 (a-h,o-z)
1613 include 'DIMENSIONS'
1614 include 'COMMON.IOUNITS'
1615 include 'COMMON.GEO'
1616 include 'COMMON.VAR'
1617 include 'COMMON.INTERACT'
1618 include 'COMMON.LOCAL'
1619 include 'COMMON.NAMES'
1620 include 'COMMON.CHAIN'
1621 include 'COMMON.FFIELD'
1622 include 'COMMON.SBRIDGE'
1623 include 'COMMON.HEADER'
1624 include 'COMMON.CONTROL'
1625 include 'COMMON.DBASE'
1626 include 'COMMON.THREAD'
1627 include 'COMMON.TIME1'
1628 C Set up variable list.
1634 if (itype(i).ne.10) then
1636 ialph(i,1)=nvar+nside
1640 if (indphi.gt.0) then
1642 else if (indback.gt.0) then
1647 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1650 c----------------------------------------------------------------------------
1651 subroutine gen_dist_constr
1652 C Generate CA distance constraints.
1653 implicit real*8 (a-h,o-z)
1654 include 'DIMENSIONS'
1655 include 'COMMON.IOUNITS'
1656 include 'COMMON.GEO'
1657 include 'COMMON.VAR'
1658 include 'COMMON.INTERACT'
1659 include 'COMMON.LOCAL'
1660 include 'COMMON.NAMES'
1661 include 'COMMON.CHAIN'
1662 include 'COMMON.FFIELD'
1663 include 'COMMON.SBRIDGE'
1664 include 'COMMON.HEADER'
1665 include 'COMMON.CONTROL'
1666 include 'COMMON.DBASE'
1667 include 'COMMON.THREAD'
1668 include 'COMMON.TIME1'
1669 dimension itype_pdb(maxres)
1670 common /pizda/ itype_pdb
1672 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1673 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1674 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1676 do i=nstart_sup,nstart_sup+nsup-1
1677 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1678 cd & ' seq_pdb', restyp(itype_pdb(i))
1679 do j=i+2,nstart_sup+nsup-1
1681 ihpb(nhpb)=i+nstart_seq-nstart_sup
1682 jhpb(nhpb)=j+nstart_seq-nstart_sup
1684 dhpb(nhpb)=dist(i,j)
1687 cd write (iout,'(a)') 'Distance constraints:'
1692 cd if (ii.gt.nres) then
1697 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1698 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1699 cd & dhpb(i),forcon(i)
1703 c----------------------------------------------------------------------------
1705 implicit real*8 (a-h,o-z)
1706 include 'DIMENSIONS'
1707 include 'COMMON.MAP'
1708 include 'COMMON.IOUNITS'
1709 character*3 angid(4) /'THE','PHI','ALP','OME'/
1710 character*80 mapcard,ucase
1712 read (inp,'(a)') mapcard
1713 mapcard=ucase(mapcard)
1714 if (index(mapcard,'PHI').gt.0) then
1716 else if (index(mapcard,'THE').gt.0) then
1718 else if (index(mapcard,'ALP').gt.0) then
1720 else if (index(mapcard,'OME').gt.0) then
1723 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1724 stop 'Error - illegal variable spec in MAP card.'
1726 call readi (mapcard,'RES1',res1(imap),0)
1727 call readi (mapcard,'RES2',res2(imap),0)
1728 if (res1(imap).eq.0) then
1729 res1(imap)=res2(imap)
1730 else if (res2(imap).eq.0) then
1731 res2(imap)=res1(imap)
1733 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1735 & 'Error - illegal definition of variable group in MAP.'
1736 stop 'Error - illegal definition of variable group in MAP.'
1738 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1739 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1740 call readi(mapcard,'NSTEP',nstep(imap),0)
1741 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1743 & 'Illegal boundary and/or step size specification in MAP.'
1744 stop 'Illegal boundary and/or step size specification in MAP.'
1749 c----------------------------------------------------------------------------
1750 csa subroutine csaread
1751 csa implicit real*8 (a-h,o-z)
1752 csa include 'DIMENSIONS'
1753 csa include 'COMMON.IOUNITS'
1754 csa include 'COMMON.GEO'
1755 csa include 'COMMON.CSA'
1756 csa include 'COMMON.BANK'
1757 csa include 'COMMON.CONTROL'
1758 csa character*80 ucase
1759 csa character*620 mcmcard
1760 csa call card_concat(mcmcard)
1762 csa call readi(mcmcard,'NCONF',nconf,50)
1763 csa call readi(mcmcard,'NADD',nadd,0)
1764 csa call readi(mcmcard,'JSTART',jstart,1)
1765 csa call readi(mcmcard,'JEND',jend,1)
1766 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1767 csa call readi(mcmcard,'N0',n0,1)
1768 csa call readi(mcmcard,'N1',n1,6)
1769 csa call readi(mcmcard,'N2',n2,4)
1770 csa call readi(mcmcard,'N3',n3,0)
1771 csa call readi(mcmcard,'N4',n4,0)
1772 csa call readi(mcmcard,'N5',n5,0)
1773 csa call readi(mcmcard,'N6',n6,10)
1774 csa call readi(mcmcard,'N7',n7,0)
1775 csa call readi(mcmcard,'N8',n8,0)
1776 csa call readi(mcmcard,'N9',n9,0)
1777 csa call readi(mcmcard,'N14',n14,0)
1778 csa call readi(mcmcard,'N15',n15,0)
1779 csa call readi(mcmcard,'N16',n16,0)
1780 csa call readi(mcmcard,'N17',n17,0)
1781 csa call readi(mcmcard,'N18',n18,0)
1783 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1785 csa call readi(mcmcard,'NDIFF',ndiff,2)
1786 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1787 csa call readi(mcmcard,'IS1',is1,1)
1788 csa call readi(mcmcard,'IS2',is2,8)
1789 csa call readi(mcmcard,'NRAN0',nran0,4)
1790 csa call readi(mcmcard,'NRAN1',nran1,2)
1791 csa call readi(mcmcard,'IRR',irr,1)
1792 csa call readi(mcmcard,'NSEED',nseed,20)
1793 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1794 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1795 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1796 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1797 csa call readi(mcmcard,'ICMAX',icmax,3)
1798 csa call readi(mcmcard,'IRESTART',irestart,0)
1799 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1802 csa call reada(mcmcard,'DELE',dele,20.0d0)
1803 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1804 csa call readi(mcmcard,'IREF',iref,0)
1805 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1806 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1807 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1808 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1809 csa write (iout,*) "NCONF_IN",nconf_in
1812 c----------------------------------------------------------------------------
1813 cfmc subroutine mcmfread
1814 cfmc implicit real*8 (a-h,o-z)
1815 cfmc include 'DIMENSIONS'
1816 cfmc include 'COMMON.MCMF'
1817 cfmc include 'COMMON.IOUNITS'
1818 cfmc include 'COMMON.GEO'
1819 cfmc character*80 ucase
1820 cfmc character*620 mcmcard
1821 cfmc call card_concat(mcmcard)
1823 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1824 cfmc write(iout,*)'MAXRANT=',maxrant
1825 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1826 cfmc write(iout,*)'MAXFAM=',maxfam
1827 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1828 cfmc write(iout,*)'NNET1=',nnet1
1829 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1830 cfmc write(iout,*)'NNET2=',nnet2
1831 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1832 cfmc write(iout,*)'NNET3=',nnet3
1833 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1834 cfmc write(iout,*)'ILASTT=',ilastt
1835 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1836 cfmc write(iout,*)'MAXSTR=',maxstr
1837 cfmc maxstr_f=maxstr/maxfam
1838 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1839 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1840 cfmc write(iout,*)'NMCMF=',nmcmf
1841 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1842 cfmc write(iout,*)'IFOCUS=',ifocus
1843 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1844 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1845 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1846 cfmc write(iout,*)'INTPRT=',intprt
1847 cfmc call readi(mcmcard,'IPRT',iprt,100)
1848 cfmc write(iout,*)'IPRT=',iprt
1849 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1850 cfmc write(iout,*)'IMAXTR=',imaxtr
1851 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1852 cfmc write(iout,*)'MAXEVEN=',maxeven
1853 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1854 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1855 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1856 cfmc write(iout,*)'INIMIN=',inimin
1857 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1858 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1859 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1860 cfmc write(iout,*)'NTHREAD=',nthread
1861 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1862 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1863 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1864 cfmc write(iout,*)'MAXPERT=',maxpert
1865 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1866 cfmc write(iout,*)'IRMSD=',irmsd
1867 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1868 cfmc write(iout,*)'DENEMIN=',denemin
1869 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1870 cfmc write(iout,*)'RCUT1S=',rcut1s
1871 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1872 cfmc write(iout,*)'RCUT1E=',rcut1e
1873 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1874 cfmc write(iout,*)'RCUT2S=',rcut2s
1875 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1876 cfmc write(iout,*)'RCUT2E=',rcut2e
1877 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1878 cfmc write(iout,*)'DPERT1=',d_pert1
1879 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1880 cfmc write(iout,*)'DPERT1A=',d_pert1a
1881 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1882 cfmc write(iout,*)'DPERT2=',d_pert2
1883 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1884 cfmc write(iout,*)'DPERT2A=',d_pert2a
1885 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1886 cfmc write(iout,*)'DPERT2B=',d_pert2b
1887 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1888 cfmc write(iout,*)'DPERT2C=',d_pert2c
1889 cfmc d_pert1=deg2rad*d_pert1
1890 cfmc d_pert1a=deg2rad*d_pert1a
1891 cfmc d_pert2=deg2rad*d_pert2
1892 cfmc d_pert2a=deg2rad*d_pert2a
1893 cfmc d_pert2b=deg2rad*d_pert2b
1894 cfmc d_pert2c=deg2rad*d_pert2c
1895 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1896 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1897 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1898 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1899 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1900 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1901 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1902 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1903 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1904 cfmc write(iout,*)'RCUTINI=',rcutini
1905 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1906 cfmc write(iout,*)'GRAT=',grat
1907 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1908 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1912 c----------------------------------------------------------------------------
1914 implicit real*8 (a-h,o-z)
1915 include 'DIMENSIONS'
1916 include 'COMMON.MCM'
1917 include 'COMMON.MCE'
1918 include 'COMMON.IOUNITS'
1920 character*320 mcmcard
1921 call card_concat(mcmcard)
1922 call readi(mcmcard,'MAXACC',maxacc,100)
1923 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1924 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1925 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1926 call readi(mcmcard,'MAXREPM',maxrepm,200)
1927 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1928 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1929 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1930 call reada(mcmcard,'E_UP',e_up,5.0D0)
1931 call reada(mcmcard,'DELTE',delte,0.1D0)
1932 call readi(mcmcard,'NSWEEP',nsweep,5)
1933 call readi(mcmcard,'NSTEPH',nsteph,0)
1934 call readi(mcmcard,'NSTEPC',nstepc,0)
1935 call reada(mcmcard,'TMIN',tmin,298.0D0)
1936 call reada(mcmcard,'TMAX',tmax,298.0D0)
1937 call readi(mcmcard,'NWINDOW',nwindow,0)
1938 call readi(mcmcard,'PRINT_MC',print_mc,0)
1939 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1940 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1941 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1942 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1943 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1944 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1945 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1946 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1947 if (nwindow.gt.0) then
1948 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1950 winlen(i)=winend(i)-winstart(i)+1
1953 if (tmax.lt.tmin) tmax=tmin
1954 if (tmax.eq.tmin) then
1958 if (nstepc.gt.0 .and. nsteph.gt.0) then
1959 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1960 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1962 C Probabilities of different move types
1963 sumpro_type(0)=0.0D0
1964 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1965 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1966 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1967 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1968 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1969 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1970 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1972 print *,'i',i,' sumprotype',sumpro_type(i)
1973 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1974 print *,'i',i,' sumprotype',sumpro_type(i)
1978 c----------------------------------------------------------------------------
1979 subroutine read_minim
1980 implicit real*8 (a-h,o-z)
1981 include 'DIMENSIONS'
1982 include 'COMMON.MINIM'
1983 include 'COMMON.IOUNITS'
1985 character*320 minimcard
1986 call card_concat(minimcard)
1987 call readi(minimcard,'MAXMIN',maxmin,2000)
1988 call readi(minimcard,'MAXFUN',maxfun,5000)
1989 call readi(minimcard,'MINMIN',minmin,maxmin)
1990 call readi(minimcard,'MINFUN',minfun,maxmin)
1991 call reada(minimcard,'TOLF',tolf,1.0D-2)
1992 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1993 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1994 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1995 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1996 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1997 & 'Options in energy minimization:'
1998 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1999 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
2000 & 'MinMin:',MinMin,' MinFun:',MinFun,
2001 & ' TolF:',TolF,' RTolF:',RTolF
2004 c----------------------------------------------------------------------------
2005 subroutine read_angles(kanal,*)
2006 implicit real*8 (a-h,o-z)
2007 include 'DIMENSIONS'
2008 include 'COMMON.GEO'
2009 include 'COMMON.VAR'
2010 include 'COMMON.CHAIN'
2011 include 'COMMON.IOUNITS'
2012 include 'COMMON.CONTROL'
2013 c Read angles from input
2015 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
2016 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
2017 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
2018 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
2021 c 9/7/01 avoid 180 deg valence angle
2022 if (theta(i).gt.179.99d0) theta(i)=179.99d0
2024 theta(i)=deg2rad*theta(i)
2025 phi(i)=deg2rad*phi(i)
2026 alph(i)=deg2rad*alph(i)
2027 omeg(i)=deg2rad*omeg(i)
2032 c----------------------------------------------------------------------------
2033 subroutine reada(rekord,lancuch,wartosc,default)
2035 character*(*) rekord,lancuch
2036 double precision wartosc,default
2039 iread=index(rekord,lancuch)
2040 if (iread.eq.0) then
2044 iread=iread+ilen(lancuch)+1
2045 read (rekord(iread:),*,err=10,end=10) wartosc
2050 c----------------------------------------------------------------------------
2051 subroutine readi(rekord,lancuch,wartosc,default)
2053 character*(*) rekord,lancuch
2054 integer wartosc,default
2057 iread=index(rekord,lancuch)
2058 if (iread.eq.0) then
2062 iread=iread+ilen(lancuch)+1
2063 read (rekord(iread:),*,err=10,end=10) wartosc
2068 c----------------------------------------------------------------------------
2069 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2072 integer tablica(dim),default
2073 character*(*) rekord,lancuch
2080 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2081 if (iread.eq.0) return
2082 iread=iread+ilen(lancuch)+1
2083 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2086 c----------------------------------------------------------------------------
2087 subroutine multreada(rekord,lancuch,tablica,dim,default)
2090 double precision tablica(dim),default
2091 character*(*) rekord,lancuch
2098 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2099 if (iread.eq.0) return
2100 iread=iread+ilen(lancuch)+1
2101 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2104 c----------------------------------------------------------------------------
2105 subroutine openunits
2106 implicit real*8 (a-h,o-z)
2107 include 'DIMENSIONS'
2110 character*16 form,nodename
2113 include 'COMMON.SETUP'
2114 include 'COMMON.IOUNITS'
2116 include 'COMMON.CONTROL'
2117 integer lenpre,lenpot,ilen,lentmp
2119 character*3 out1file_text,ucase
2122 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2123 call getenv_loc("PREFIX",prefix)
2125 call getenv_loc("POT",pot)
2126 call getenv_loc("DIRTMP",tmpdir)
2127 call getenv_loc("CURDIR",curdir)
2128 call getenv_loc("OUT1FILE",out1file_text)
2129 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2130 out1file_text=ucase(out1file_text)
2131 if (out1file_text(1:1).eq."Y") then
2134 out1file=fg_rank.gt.0
2139 if (lentmp.gt.0) then
2140 write (*,'(80(1h!))')
2141 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2142 write (*,'(80(1h!))')
2143 write (*,*)"All output files will be on node /tmp directory."
2145 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2146 if (me.eq.king) then
2147 write (*,*) "The master node is ",nodename
2148 else if (fg_rank.eq.0) then
2149 write (*,*) "I am the CG slave node ",nodename
2151 write (*,*) "I am the FG slave node ",nodename
2154 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2155 lenpre = lentmp+lenpre+1
2157 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2158 C Get the names and open the input files
2159 #if defined(WINIFL) || defined(WINPGI)
2160 open(1,file=pref_orig(:ilen(pref_orig))//
2161 & '.inp',status='old',readonly,shared)
2162 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2163 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2164 C Get parameter filenames and open the parameter files.
2165 call getenv_loc('BONDPAR',bondname)
2166 open (ibond,file=bondname,status='old',readonly,shared)
2167 call getenv_loc('THETPAR',thetname)
2168 open (ithep,file=thetname,status='old',readonly,shared)
2170 call getenv_loc('THETPARPDB',thetname_pdb)
2171 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2173 call getenv_loc('ROTPAR',rotname)
2174 open (irotam,file=rotname,status='old',readonly,shared)
2176 call getenv_loc('ROTPARPDB',rotname_pdb)
2177 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2179 call getenv_loc('TORPAR',torname)
2180 open (itorp,file=torname,status='old',readonly,shared)
2181 call getenv_loc('TORDPAR',tordname)
2182 open (itordp,file=tordname,status='old',readonly,shared)
2183 call getenv_loc('FOURIER',fouriername)
2184 open (ifourier,file=fouriername,status='old',readonly,shared)
2185 call getenv_loc('ELEPAR',elename)
2186 open (ielep,file=elename,status='old',readonly,shared)
2187 call getenv_loc('SIDEPAR',sidename)
2188 open (isidep,file=sidename,status='old',readonly,shared)
2189 #elif (defined CRAY) || (defined AIX)
2190 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2192 c print *,"Processor",myrank," opened file 1"
2193 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2194 c print *,"Processor",myrank," opened file 9"
2195 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2196 C Get parameter filenames and open the parameter files.
2197 call getenv_loc('BONDPAR',bondname)
2198 open (ibond,file=bondname,status='old',action='read')
2199 c print *,"Processor",myrank," opened file IBOND"
2200 call getenv_loc('THETPAR',thetname)
2201 open (ithep,file=thetname,status='old',action='read')
2202 c print *,"Processor",myrank," opened file ITHEP"
2204 call getenv_loc('THETPARPDB',thetname_pdb)
2205 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2207 call getenv_loc('ROTPAR',rotname)
2208 open (irotam,file=rotname,status='old',action='read')
2209 c print *,"Processor",myrank," opened file IROTAM"
2211 call getenv_loc('ROTPARPDB',rotname_pdb)
2212 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2214 call getenv_loc('TORPAR',torname)
2215 open (itorp,file=torname,status='old',action='read')
2216 c print *,"Processor",myrank," opened file ITORP"
2217 call getenv_loc('TORDPAR',tordname)
2218 open (itordp,file=tordname,status='old',action='read')
2219 c print *,"Processor",myrank," opened file ITORDP"
2220 call getenv_loc('SCCORPAR',sccorname)
2221 open (isccor,file=sccorname,status='old',action='read')
2222 c print *,"Processor",myrank," opened file ISCCOR"
2223 call getenv_loc('FOURIER',fouriername)
2224 open (ifourier,file=fouriername,status='old',action='read')
2225 c print *,"Processor",myrank," opened file IFOURIER"
2226 call getenv_loc('ELEPAR',elename)
2227 open (ielep,file=elename,status='old',action='read')
2228 c print *,"Processor",myrank," opened file IELEP"
2229 call getenv_loc('SIDEPAR',sidename)
2230 open (isidep,file=sidename,status='old',action='read')
2231 c print *,"Processor",myrank," opened file ISIDEP"
2232 c print *,"Processor",myrank," opened parameter files"
2234 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2235 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2236 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2237 C Get parameter filenames and open the parameter files.
2238 call getenv_loc('BONDPAR',bondname)
2239 open (ibond,file=bondname,status='old')
2240 call getenv_loc('THETPAR',thetname)
2241 open (ithep,file=thetname,status='old')
2243 call getenv_loc('THETPARPDB',thetname_pdb)
2244 open (ithep_pdb,file=thetname_pdb,status='old')
2246 call getenv_loc('ROTPAR',rotname)
2247 open (irotam,file=rotname,status='old')
2249 call getenv_loc('ROTPARPDB',rotname_pdb)
2250 open (irotam_pdb,file=rotname_pdb,status='old')
2252 call getenv_loc('TORPAR',torname)
2253 open (itorp,file=torname,status='old')
2254 call getenv_loc('TORDPAR',tordname)
2255 open (itordp,file=tordname,status='old')
2256 call getenv_loc('SCCORPAR',sccorname)
2257 open (isccor,file=sccorname,status='old')
2258 call getenv_loc('FOURIER',fouriername)
2259 open (ifourier,file=fouriername,status='old')
2260 call getenv_loc('ELEPAR',elename)
2261 open (ielep,file=elename,status='old')
2262 call getenv_loc('SIDEPAR',sidename)
2263 open (isidep,file=sidename,status='old')
2265 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2267 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2268 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2269 C Get parameter filenames and open the parameter files.
2270 call getenv_loc('BONDPAR',bondname)
2271 open (ibond,file=bondname,status='old',action='read')
2272 call getenv_loc('THETPAR',thetname)
2273 open (ithep,file=thetname,status='old',action='read')
2275 call getenv_loc('THETPARPDB',thetname_pdb)
2276 print *,"thetname_pdb ",thetname_pdb
2277 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2278 print *,ithep_pdb," opened"
2280 call getenv_loc('ROTPAR',rotname)
2281 open (irotam,file=rotname,status='old',action='read')
2283 call getenv_loc('ROTPARPDB',rotname_pdb)
2284 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2286 call getenv_loc('TORPAR',torname)
2287 open (itorp,file=torname,status='old',action='read')
2288 call getenv_loc('TORDPAR',tordname)
2289 open (itordp,file=tordname,status='old',action='read')
2290 call getenv_loc('SCCORPAR',sccorname)
2291 open (isccor,file=sccorname,status='old',action='read')
2292 call getenv_loc('FOURIER',fouriername)
2293 open (ifourier,file=fouriername,status='old',action='read')
2294 call getenv_loc('ELEPAR',elename)
2295 open (ielep,file=elename,status='old',action='read')
2296 call getenv_loc('SIDEPAR',sidename)
2297 open (isidep,file=sidename,status='old',action='read')
2301 C 8/9/01 In the newest version SCp interaction constants are read from a file
2302 C Use -DOLDSCP to use hard-coded constants instead.
2304 call getenv_loc('SCPPAR',scpname)
2305 #if defined(WINIFL) || defined(WINPGI)
2306 open (iscpp,file=scpname,status='old',readonly,shared)
2307 #elif (defined CRAY) || (defined AIX)
2308 open (iscpp,file=scpname,status='old',action='read')
2310 open (iscpp,file=scpname,status='old')
2312 open (iscpp,file=scpname,status='old',action='read')
2315 call getenv_loc('PATTERN',patname)
2316 #if defined(WINIFL) || defined(WINPGI)
2317 open (icbase,file=patname,status='old',readonly,shared)
2318 #elif (defined CRAY) || (defined AIX)
2319 open (icbase,file=patname,status='old',action='read')
2321 open (icbase,file=patname,status='old')
2323 open (icbase,file=patname,status='old',action='read')
2326 C Open output file only for CG processes
2327 c print *,"Processor",myrank," fg_rank",fg_rank
2328 if (fg_rank.eq.0) then
2330 if (nodes.eq.1) then
2333 npos = dlog10(dfloat(nodes-1))+1
2335 if (npos.lt.3) npos=3
2336 write (liczba,'(i1)') npos
2337 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2339 write (liczba,form) me
2340 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2341 & liczba(:ilen(liczba))
2342 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2344 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2346 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2347 & liczba(:ilen(liczba))//'.mol2'
2348 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//
2349 & liczba(:ilen(liczba))//'.x'
2350 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2351 & liczba(:ilen(liczba))//'.stat'
2353 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2354 & //liczba(:ilen(liczba))//'.stat')
2355 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2358 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2359 & liczba(:ilen(liczba))//'.const'
2364 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2365 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2366 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2367 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2368 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x'
2369 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2371 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2373 rest2name=prefix(:ilen(prefix))//'.rst'
2375 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2378 #if defined(AIX) || defined(PGI)
2379 if (me.eq.king .or. .not. out1file)
2380 & open(iout,file=outname,status='unknown')
2383 if (fg_rank.gt.0) then
2384 write (liczba,'(i3.3)') myrank/nfgtasks
2385 write (ll,'(bz,i3.3)') fg_rank
2386 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2392 open(igeom,file=intname,status='unknown',position='append')
2393 open(ipdb,file=pdbname,status='unknown')
2394 open(imol2,file=mol2name,status='unknown')
2395 open(istat,file=statname,status='unknown',position='append')
2397 c1out open(iout,file=outname,status='unknown')
2400 if (me.eq.king .or. .not.out1file)
2401 & open(iout,file=outname,status='unknown')
2404 if (fg_rank.gt.0) then
2405 print "Processor",fg_rank," opening output file"
2406 write (liczba,'(i3.3)') myrank/nfgtasks
2407 write (ll,'(bz,i3.3)') fg_rank
2408 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2414 open(igeom,file=intname,status='unknown',access='append')
2415 open(ipdb,file=pdbname,status='unknown')
2416 open(imol2,file=mol2name,status='unknown')
2417 open(istat,file=statname,status='unknown',access='append')
2419 c1out open(iout,file=outname,status='unknown')
2422 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2423 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2424 csa csa_history=prefix(:lenpre)//'.CSA.history'
2425 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2426 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2427 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2428 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2429 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2430 csa csa_int=prefix(:lenpre)//'.int'
2431 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2432 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2433 csa csa_in=prefix(:lenpre)//'.CSA.in'
2434 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2437 write (iout,'(80(1h-))')
2438 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2439 write (iout,'(80(1h-))')
2440 write (iout,*) "Input file : ",
2441 & pref_orig(:ilen(pref_orig))//'.inp'
2442 write (iout,*) "Output file : ",
2443 & outname(:ilen(outname))
2445 write (iout,*) "Sidechain potential file : ",
2446 & sidename(:ilen(sidename))
2448 write (iout,*) "SCp potential file : ",
2449 & scpname(:ilen(scpname))
2451 write (iout,*) "Electrostatic potential file : ",
2452 & elename(:ilen(elename))
2453 write (iout,*) "Cumulant coefficient file : ",
2454 & fouriername(:ilen(fouriername))
2455 write (iout,*) "Torsional parameter file : ",
2456 & torname(:ilen(torname))
2457 write (iout,*) "Double torsional parameter file : ",
2458 & tordname(:ilen(tordname))
2459 write (iout,*) "SCCOR parameter file : ",
2460 & sccorname(:ilen(sccorname))
2461 write (iout,*) "Bond & inertia constant file : ",
2462 & bondname(:ilen(bondname))
2463 write (iout,*) "Bending parameter file : ",
2464 & thetname(:ilen(thetname))
2465 write (iout,*) "Rotamer parameter file : ",
2466 & rotname(:ilen(rotname))
2467 write (iout,*) "Threading database : ",
2468 & patname(:ilen(patname))
2470 &write (iout,*)" DIRTMP : ",
2472 write (iout,'(80(1h-))')
2476 c----------------------------------------------------------------------------
2477 subroutine card_concat(card)
2478 implicit real*8 (a-h,o-z)
2479 include 'DIMENSIONS'
2480 include 'COMMON.IOUNITS'
2482 character*80 karta,ucase
2484 read (inp,'(a)') karta
2487 do while (karta(80:80).eq.'&')
2488 card=card(:ilen(card)+1)//karta(:79)
2489 read (inp,'(a)') karta
2492 card=card(:ilen(card)+1)//karta
2495 c----------------------------------------------------------------------------------
2497 implicit real*8 (a-h,o-z)
2498 include 'DIMENSIONS'
2499 include 'COMMON.CHAIN'
2500 include 'COMMON.IOUNITS'
2502 include 'COMMON.CONTROL'
2503 open(irest2,file=rest2name,status='unknown')
2504 read(irest2,*) totT,EK,potE,totE,t_bath
2506 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2509 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2511 if(usampl.or.homol_nset.gt.1) then
2512 read (irest2,*) iset
2517 c---------------------------------------------------------------------------------
2518 subroutine read_fragments
2519 implicit real*8 (a-h,o-z)
2520 include 'DIMENSIONS'
2524 include 'COMMON.SETUP'
2525 include 'COMMON.CHAIN'
2526 include 'COMMON.IOUNITS'
2528 include 'COMMON.CONTROL'
2530 read(inp,*) nset,nfrag,npair,nfrag_back
2531 if(me.eq.king.or..not.out1file)
2532 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2533 & " nfrag_back",nfrag_back
2535 read(inp,*) mset(iset1)
2537 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2539 if(me.eq.king.or..not.out1file)
2540 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2541 & ifrag(2,i,iset1), qinfrag(i,iset1)
2544 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2546 if(me.eq.king.or..not.out1file)
2547 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2548 & ipair(2,i,iset1), qinpair(i,iset1)
2551 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2552 & wfrag_back(3,i,iset1),
2553 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2554 if(me.eq.king.or..not.out1file)
2555 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2556 & wfrag_back(2,i,iset1),
2557 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2558 & ifrag_back(2,i,iset1)
2563 c-------------------------------------------------------------------------------
2564 subroutine read_dist_constr
2565 implicit real*8 (a-h,o-z)
2566 include 'DIMENSIONS'
2570 include 'COMMON.SETUP'
2571 include 'COMMON.CONTROL'
2572 include 'COMMON.CHAIN'
2573 include 'COMMON.IOUNITS'
2574 include 'COMMON.SBRIDGE'
2575 integer ifrag_(2,100),ipair_(2,100)
2576 double precision wfrag_(100),wpair_(100)
2577 character*500 controlcard
2578 c write (iout,*) "Calling read_dist_constr"
2579 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2581 call card_concat(controlcard)
2582 call readi(controlcard,"NFRAG",nfrag_,0)
2583 call readi(controlcard,"NPAIR",npair_,0)
2584 call readi(controlcard,"NDIST",ndist_,0)
2585 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2587 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2588 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2589 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2590 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2591 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2592 c write (iout,*) "IFRAG"
2594 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2596 c write (iout,*) "IPAIR"
2598 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2600 if (.not.refstr .and. nfrag.gt.0) then
2602 & "ERROR: no reference structure to compute distance restraints"
2604 & "Restraints must be specified explicitly (NDIST=number)"
2607 if (nfrag.lt.2 .and. npair.gt.0) then
2608 write (iout,*) "ERROR: Less than 2 fragments specified",
2609 & " but distance restraints between pairs requested"
2614 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2615 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2616 & ifrag_(2,i)=nstart_sup+nsup-1
2617 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2619 if (wfrag_(i).gt.0.0d0) then
2620 do j=ifrag_(1,i),ifrag_(2,i)-1
2621 do k=j+1,ifrag_(2,i)
2622 c write (iout,*) "j",j," k",k
2624 if (constr_dist.eq.1) then
2629 forcon(nhpb)=wfrag_(i)
2630 else if (constr_dist.eq.2) then
2631 if (ddjk.le.dist_cut) then
2636 forcon(nhpb)=wfrag_(i)
2643 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2646 if (.not.out1file .or. me.eq.king)
2647 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2648 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2650 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2651 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2658 if (wpair_(i).gt.0.0d0) then
2666 do j=ifrag_(1,ii),ifrag_(2,ii)
2667 do k=ifrag_(1,jj),ifrag_(2,jj)
2671 forcon(nhpb)=wpair_(i)
2672 dhpb(nhpb)=dist(j,k)
2674 if (.not.out1file .or. me.eq.king)
2675 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2676 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2678 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2679 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2686 if (constr_dist.eq.11) then
2687 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2688 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2689 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2691 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2692 & ibecarb(i),forcon(nhpb+1)
2694 if (forcon(nhpb+1).gt.0.0d0) then
2696 if (ibecarb(i).gt.0) then
2697 ihpb(i)=ihpb(i)+nres
2698 jhpb(i)=jhpb(i)+nres
2700 if (dhpb(nhpb).eq.0.0d0)
2701 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2705 if (.not.out1file .or. me.eq.king) then
2708 if (constr_dist.eq.11) then
2709 write (iout,'(a,3i5,2f8.2,i2,2f10.1)') "+dist.constr11 ",
2710 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i),
2713 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2714 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2723 c-------------------------------------------------------------------------------
2725 subroutine read_constr_homology
2727 include 'DIMENSIONS'
2731 include 'COMMON.SETUP'
2732 include 'COMMON.CONTROL'
2733 include 'COMMON.CHAIN'
2734 include 'COMMON.IOUNITS'
2736 include 'COMMON.GEO'
2737 include 'COMMON.INTERACT'
2738 include 'COMMON.NAMES'
2740 c For new homol impl
2742 include 'COMMON.VAR'
2745 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2747 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2748 c & sigma_odl_temp(maxres,maxres,max_template)
2750 character*24 model_ki_dist, model_ki_angle
2751 character*500 controlcard
2752 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2753 logical lprn /.true./
2758 c FP - Nov. 2014 Temporary specifications for new vars
2760 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
2761 double precision, dimension (max_template,maxres) :: rescore
2762 double precision, dimension (max_template,maxres) :: rescore2
2763 character*24 pdbfile,tpl_k_rescore
2764 c -----------------------------------------------------------------
2765 c Reading multiple PDB ref structures and calculation of retraints
2766 c not using pre-computed ones stored in files model_ki_{dist,angle}
2768 c -----------------------------------------------------------------
2771 c Alternative: reading from input
2772 call card_concat(controlcard)
2773 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2774 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2775 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2776 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2777 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2778 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2779 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2780 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2781 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2782 if(.not.read2sigma.and.start_from_model) then
2783 write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2784 start_from_model=.false.
2786 if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON'
2787 if(start_from_model .and. rest) then
2788 write(iout,*) 'START_FROM_MODELS is OFF'
2789 write(iout,*) 'remove restart keyword from input'
2791 if (homol_nset.gt.1)then
2792 call card_concat(controlcard)
2793 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2794 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2795 write(iout,*) "iset homology_weight "
2797 write(iout,*) i,waga_homology(i)
2800 iset=mod(kolor,homol_nset)+1
2803 waga_homology(1)=1.0
2806 cd write (iout,*) "nnt",nnt," nct",nct
2813 write(iout,*) 'nnt=',nnt,'nct=',nct
2816 do k=1,constr_homology
2829 if (read_homol_frag) then
2830 call read_klapaucjusz
2833 do k=1,constr_homology
2835 read(inp,'(a)') pdbfile
2836 c Next stament causes error upon compilation (?)
2837 c if(me.eq.king.or. .not. out1file)
2838 c write (iout,'(2a)') 'PDB data will be read from file ',
2839 c & pdbfile(:ilen(pdbfile))
2840 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2841 & pdbfile(:ilen(pdbfile))
2842 open(ipdbin,file=pdbfile,status='old',err=33)
2844 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2845 & pdbfile(:ilen(pdbfile))
2848 c print *,'Begin reading pdb data'
2850 c Files containing res sim or local scores (former containing sigmas)
2853 write(kic2,'(bz,i2.2)') k
2855 tpl_k_rescore="template"//kic2//".sco"
2858 if (read2sigma) then
2859 call readpdb_template(k)
2864 c Distance restraints
2867 C Copy the coordinates from reference coordinates (?)
2871 c write (iout,*) "c(",j,i,") =",c(j,i)
2875 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2877 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2878 open (ientin,file=tpl_k_rescore,status='old')
2879 if (nnt.gt.1) rescore(k,1)=0.0d0
2880 do irec=nnt,nct ! loop for reading res sim
2881 if (read2sigma) then
2882 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2885 idomain(k,i_tmp)=idomain_tmp
2886 rescore(k,i_tmp)=rescore_tmp
2887 rescore2(k,i_tmp)=rescore2_tmp
2888 write(iout,'(a7,i5,2f10.5,i5)') "rescore",
2889 & i_tmp,rescore2_tmp,rescore_tmp,
2893 read (ientin,*,end=1401) rescore_tmp
2895 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2896 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2897 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2902 if (waga_dist.ne.0.0d0) then
2910 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2911 c write (iout,*) k,i,j,distal,dist2_cut
2913 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2914 & .and. distal.le.dist2_cut ) then
2920 c write (iout,*) "k",k
2921 c write (iout,*) "i",i," j",j," constr_homology",
2926 if (read2sigma) then
2929 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2931 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2932 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2933 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2935 if (odl(k,ii).le.dist_cut) then
2936 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2939 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2940 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2942 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2943 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2947 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2950 l_homo(k,ii)=.false.
2957 c Theta, dihedral and SC retraints
2959 if (waga_angle.gt.0.0d0) then
2960 c open (ientin,file=tpl_k_sigma_dih,status='old')
2961 c do irec=1,maxres-3 ! loop for reading sigma_dih
2962 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2963 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2964 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2965 c & sigma_dih(k,i+nnt-1)
2970 if (idomain(k,i).eq.0) then
2974 dih(k,i)=phiref(i) ! right?
2975 c read (ientin,*) sigma_dih(k,i) ! original variant
2976 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2977 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2978 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2979 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2980 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2982 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2983 & rescore(k,i-2)+rescore(k,i-3))/4.0
2984 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2985 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2986 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2987 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2988 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2989 c Instead of res sim other local measure of b/b str reliability possible
2990 if (sigma_dih(k,i).ne.0)
2991 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2992 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2997 if (waga_theta.gt.0.0d0) then
2998 c open (ientin,file=tpl_k_sigma_theta,status='old')
2999 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3000 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3001 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3002 c & sigma_theta(k,i+nnt-1)
3007 do i = nnt+2,nct ! right? without parallel.
3008 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3009 c do i=ithet_start,ithet_end ! with FG parallel.
3010 if (idomain(k,i).eq.0) then
3011 sigma_theta(k,i)=0.0
3014 thetatpl(k,i)=thetaref(i)
3015 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3016 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3017 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3018 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3019 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3020 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3021 & rescore(k,i-2))/3.0
3022 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3023 if (sigma_theta(k,i).ne.0)
3024 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3026 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3027 c rescore(k,i-2) ! right expression ?
3028 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3032 if (waga_d.gt.0.0d0) then
3033 c open (ientin,file=tpl_k_sigma_d,status='old')
3034 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3035 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3036 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3037 c & sigma_d(k,i+nnt-1)
3041 do i = nnt,nct ! right? without parallel.
3042 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3043 c do i=loc_start,loc_end ! with FG parallel.
3044 if (itype(i).eq.10) cycle
3045 if (idomain(k,i).eq.0 ) then
3052 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3053 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3054 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3055 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3056 sigma_d(k,i)=rescore(k,i) ! right expression ?
3057 if (sigma_d(k,i).ne.0)
3058 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3060 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3061 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3062 c read (ientin,*) sigma_d(k,i) ! 1st variant
3063 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
3068 c remove distance restraints not used in any model from the list
3069 c shift data in all arrays
3071 if (waga_dist.ne.0.0d0) then
3077 if (ii_in_use(ii).eq.0.and.liiflag) then
3081 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3082 & .not.liiflag.and.ii.eq.lim_odl) then
3083 if (ii.eq.lim_odl) then
3084 iishift=ii-iistart+1
3089 do ki=iistart,lim_odl-iishift
3090 ires_homo(ki)=ires_homo(ki+iishift)
3091 jres_homo(ki)=jres_homo(ki+iishift)
3092 ii_in_use(ki)=ii_in_use(ki+iishift)
3093 do k=1,constr_homology
3094 odl(k,ki)=odl(k,ki+iishift)
3095 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3096 l_homo(k,ki)=l_homo(k,ki+iishift)
3100 lim_odl=lim_odl-iishift
3106 endif ! .not. klapaucjusz
3108 if (constr_homology.gt.0) call homology_partition
3109 if (constr_homology.gt.0) call init_int_table
3110 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
3111 cd & "lim_xx=",lim_xx
3112 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3113 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3117 if (.not.lprn) return
3118 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3119 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3120 write (iout,*) "Distance restraints from templates"
3122 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3123 & ii,ires_homo(ii),jres_homo(ii),
3124 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3125 & ki=1,constr_homology)
3127 write (iout,*) "Dihedral angle restraints from templates"
3129 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3130 & (rad2deg*dih(ki,i),
3131 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3133 write (iout,*) "Virtual-bond angle restraints from templates"
3135 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3136 & (rad2deg*thetatpl(ki,i),
3137 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3139 write (iout,*) "SC restraints from templates"
3141 write(iout,'(i5,100(4f8.2,4x))') i,
3142 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3143 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3146 c -----------------------------------------------------------------
3149 c -----------------------------------------------------------------
3150 subroutine read_klapaucjusz
3152 include 'DIMENSIONS'
3156 include 'COMMON.SETUP'
3157 include 'COMMON.CONTROL'
3158 include 'COMMON.CHAIN'
3159 include 'COMMON.IOUNITS'
3161 include 'COMMON.GEO'
3162 include 'COMMON.INTERACT'
3163 include 'COMMON.NAMES'
3164 character*256 fragfile
3165 integer ninclust(maxclust),inclust(max_template,maxclust),
3166 & nresclust(maxclust),iresclust(maxres,maxclust)
3169 character*24 model_ki_dist, model_ki_angle
3170 character*500 controlcard
3171 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3172 logical lprn /.true./
3178 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3179 double precision, dimension (max_template,maxres) :: rescore
3180 double precision, dimension (max_template,maxres) :: rescore2
3181 character*24 pdbfile,tpl_k_rescore
3184 c For new homol impl
3186 include 'COMMON.VAR'
3188 call getenv("FRAGFILE",fragfile)
3189 open(ientin,file=fragfile,status="old",err=10)
3190 read(ientin,*) constr_homology,nclust
3196 do k=1,constr_homology
3197 read(ientin,'(a)') pdbfile
3198 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3199 & pdbfile(:ilen(pdbfile))
3200 open(ipdbin,file=pdbfile,status='old',err=33)
3202 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3203 & pdbfile(:ilen(pdbfile))
3207 call readpdb_template(k)
3215 read(ientin,*) ninclust(i),nresclust(i)
3216 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3217 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3220 c Loop over clusters
3223 do ll = 1,ninclust(l)
3231 idomain(k,iresclust(i,l)+1) = 1
3233 idomain(k,iresclust(i,l)) = 1
3237 c Distance restraints
3240 C Copy the coordinates from reference coordinates (?)
3244 c write (iout,*) "c(",j,i,") =",c(j,i)
3247 call int_from_cart(.true.,.false.)
3248 call sc_loc_geom(.false.)
3250 thetaref(i)=theta(i)
3253 if (waga_dist.ne.0.0d0) then
3261 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3262 c write (iout,*) k,i,j,distal,dist2_cut
3264 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3265 & .and. distal.le.dist2_cut ) then
3271 c write (iout,*) "k",k
3272 c write (iout,*) "i",i," j",j," constr_homology",
3277 if (read2sigma) then
3280 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3282 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3283 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3284 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3286 if (odl(k,ii).le.dist_cut) then
3287 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3290 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3291 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3293 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3294 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3298 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3301 c l_homo(k,ii)=.false.
3308 c Theta, dihedral and SC retraints
3310 if (waga_angle.gt.0.0d0) then
3312 if (idomain(k,i).eq.0) then
3313 c sigma_dih(k,i)=0.0
3317 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3318 & rescore(k,i-2)+rescore(k,i-3))/4.0
3319 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3320 c & " sigma_dihed",sigma_dih(k,i)
3321 if (sigma_dih(k,i).ne.0)
3322 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3327 if (waga_theta.gt.0.0d0) then
3329 if (idomain(k,i).eq.0) then
3330 c sigma_theta(k,i)=0.0
3333 thetatpl(k,i)=thetaref(i)
3334 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3335 & rescore(k,i-2))/3.0
3336 if (sigma_theta(k,i).ne.0)
3337 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3341 if (waga_d.gt.0.0d0) then
3343 if (itype(i).eq.10) cycle
3344 if (idomain(k,i).eq.0 ) then
3351 sigma_d(k,i)=rescore(k,i)
3352 if (sigma_d(k,i).ne.0)
3353 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3354 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3360 c remove distance restraints not used in any model from the list
3361 c shift data in all arrays
3363 if (waga_dist.ne.0.0d0) then
3369 if (ii_in_use(ii).eq.0.and.liiflag) then
3373 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3374 & .not.liiflag.and.ii.eq.lim_odl) then
3375 if (ii.eq.lim_odl) then
3376 iishift=ii-iistart+1
3381 do ki=iistart,lim_odl-iishift
3382 ires_homo(ki)=ires_homo(ki+iishift)
3383 jres_homo(ki)=jres_homo(ki+iishift)
3384 ii_in_use(ki)=ii_in_use(ki+iishift)
3385 do k=1,constr_homology
3386 odl(k,ki)=odl(k,ki+iishift)
3387 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3388 l_homo(k,ki)=l_homo(k,ki+iishift)
3392 lim_odl=lim_odl-iishift
3399 10 stop "Error infragment file"
3401 c----------------------------------------------------------------------
3404 subroutine flush(iu)
3409 subroutine flush(iu)
3415 c------------------------------------------------------------------------------
3416 subroutine copy_to_tmp(source)
3417 include "DIMENSIONS"
3418 include "COMMON.IOUNITS"
3419 character*(*) source
3420 character* 256 tmpfile
3424 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3425 inquire(file=tmpfile,exist=ex)
3427 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3428 & " to temporary directory..."
3429 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3430 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3434 c------------------------------------------------------------------------------
3435 subroutine move_from_tmp(source)
3436 include "DIMENSIONS"
3437 include "COMMON.IOUNITS"
3438 character*(*) source
3441 write (*,*) "Moving ",source(:ilen(source)),
3442 & " from temporary directory to working directory"
3443 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3444 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3447 c------------------------------------------------------------------------------
3448 subroutine random_init(seed)
3450 C Initialize random number generator
3452 implicit real*8 (a-h,o-z)
3453 include 'DIMENSIONS'
3459 logical OKRandom, prng_restart
3461 integer iseed_array(4)
3463 include 'COMMON.IOUNITS'
3464 include 'COMMON.TIME1'
3465 include 'COMMON.THREAD'
3466 include 'COMMON.SBRIDGE'
3467 include 'COMMON.CONTROL'
3468 include 'COMMON.MCM'
3469 include 'COMMON.MAP'
3470 include 'COMMON.HEADER'
3471 csa include 'COMMON.CSA'
3472 include 'COMMON.CHAIN'
3473 include 'COMMON.MUCA'
3475 include 'COMMON.FFIELD'
3476 include 'COMMON.SETUP'
3477 iseed=-dint(dabs(seed))
3478 if (iseed.eq.0) then
3479 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3480 & 'Random seed undefined. The program will stop.'
3481 write (*,'(/80(1h*)/20x,a/80(1h*))')
3482 & 'Random seed undefined. The program will stop.'
3484 call mpi_finalize(mpi_comm_world,ierr)
3486 stop 'Bad random seed.'
3489 if (fg_rank.eq.0) then
3493 if(me.eq.king .or. .not. out1file)
3494 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3495 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3496 OKRandom = prng_restart(me,iseedi8)
3499 tmp=65536.0d0**(4-i)
3500 iseed_array(i) = dint(seed/tmp)
3501 seed=seed-iseed_array(i)*tmp
3503 if(me.eq.king .or. .not. out1file)
3504 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3505 & (iseed_array(i),i=1,4)
3506 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3507 & (iseed_array(i),i=1,4)
3508 OKRandom = prng_restart(me,iseed_array)
3511 r1=ran_number(0.0D0,1.0D0)
3512 if(me.eq.king .or. .not. out1file)
3513 & write (iout,*) 'ran_num',r1
3514 if (r1.lt.0.0d0) OKRandom=.false.
3516 if (.not.OKRandom) then
3517 write (iout,*) 'PRNG IS NOT WORKING!!!'
3518 print *,'PRNG IS NOT WORKING!!!'
3521 call mpi_abort(mpi_comm_world,error_msg,ierr)
3524 write (iout,*) 'too many processors for parallel prng'
3525 write (*,*) 'too many processors for parallel prng'
3533 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)