8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
14 C Read job setup parameters
16 C Read force-field parameters except weights
18 C Read control parameters for energy minimzation if required
19 if (minim) call read_minim
20 C Read MCM control parameters if required
21 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
22 C Read MD control parameters if reqjuired
23 if (modecalc.eq.12) call read_MDpar
24 C Read MREMD control parameters if required
25 if (modecalc.eq.14) then
29 C Read MUCA control parameters if required
30 if (lmuca) call read_muca
31 C Read CSA control parameters if required (from fort.40 if exists
32 C otherwise from general input file)
33 if (modecalc.eq.8) then
34 inquire (file="fort.40",exist=file_exist)
35 if (.not.file_exist) call csaread
37 cfmc if (modecalc.eq.10) call mcmfread
38 c Read energy-term weights and disulfide parameters
40 C Read molecule information, molecule geometry, energy-term weights, and
41 C restraints if requested
43 C Print restraint information
45 if (.not. out1file .or. me.eq.king) then
48 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
49 & "The following",nhpb-nss,
50 & " distance restraints have been imposed:",
51 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
54 write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
55 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
60 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
61 & "The following",npeak,
62 & " NMR peak restraints have been imposed:",
63 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
64 & " score"," type"," ipeak"
66 do j=ipeak(1,i),ipeak(2,i)
67 write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
68 & jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
69 & forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
72 write (iout,*) "The ipeak array"
74 write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
80 c print *,"Processor",myrank," leaves READRTNS"
83 C-------------------------------------------------------------------------------
84 subroutine read_control
92 logical OKRandom, prng_restart
95 include 'COMMON.IOUNITS'
96 include 'COMMON.TIME1'
97 include 'COMMON.THREAD'
98 include 'COMMON.SBRIDGE'
99 include 'COMMON.CONTROL'
100 include 'COMMON.SAXS'
103 include 'COMMON.HEADER'
105 include 'COMMON.CHAIN'
106 include 'COMMON.MUCA'
108 include 'COMMON.FFIELD'
109 include 'COMMON.INTERACT'
110 include 'COMMON.SETUP'
111 include 'COMMON.SPLITELE'
112 include 'COMMON.SHIELD'
115 integer KDIAG,ICORFL,IXDR
116 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
117 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
119 character*320 controlcard
120 double precision seed
125 read (INP,'(a)') titel
126 call card_concat(controlcard)
127 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
128 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
129 call reada(controlcard,'SEED',seed,0.0D0)
130 call random_init(seed)
131 C Set up the time limit (caution! The time must be input in minutes!)
132 read_cart=index(controlcard,'READ_CART').gt.0
133 out_cart=index(controlcard,'OUT_CART').gt.0
134 out_int=index(controlcard,'OUT_INT').gt.0
135 gmatout=index(controlcard,'GMATOUT').gt.0
136 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
137 C this variable with_theta_constr is the variable which allow to read and execute the
138 C constrains on theta angles WITH_THETA_CONSTR is the keyword
139 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
140 write (iout,*) "with_theta_constr ",with_theta_constr
141 call readi(controlcard,'NSAXS',nsaxs,0)
142 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
143 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
144 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
145 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
146 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
147 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
148 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
149 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
150 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
151 call readi(controlcard,'SYM',symetr,1)
152 call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
153 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
154 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
155 c call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
156 c call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
157 c call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
158 c call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
159 c call reada(controlcard,'DRMS',drms,0.1D0)
160 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
161 c write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
162 c write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
163 c write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
164 c write (iout,'(a,f10.1)')'DRMS = ',drms
165 cc write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
166 c write (iout,'(a,f10.1)') 'Time limit (min):',timlim
168 call readi(controlcard,'NZ_START',nz_start,0)
169 call readi(controlcard,'NZ_END',nz_end,0)
170 c call readi(controlcard,'IZ_SC',iz_sc,0)
172 safety = 60.0d0*safety
174 call reada(controlcard,"T_BATH",t_bath,300.0d0)
175 minim=(index(controlcard,'MINIMIZE').gt.0)
176 dccart=(index(controlcard,'CART').gt.0)
177 overlapsc=(index(controlcard,'OVERLAP').gt.0)
178 overlapsc=.not.overlapsc
179 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
180 searchsc=.not.searchsc
181 sideadd=(index(controlcard,'SIDEADD').gt.0)
182 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
183 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
184 outpdb=(index(controlcard,'PDBOUT').gt.0)
185 outmol2=(index(controlcard,'MOL2OUT').gt.0)
186 pdbref=(index(controlcard,'PDBREF').gt.0)
187 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
188 indpdb=index(controlcard,'PDBSTART')
189 extconf=(index(controlcard,'EXTCONF').gt.0)
190 AFMlog=(index(controlcard,'AFM'))
191 selfguide=(index(controlcard,'SELFGUIDE'))
192 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
193 call readi(controlcard,'TUBEMOD',tubelog,0)
194 c write (iout,*) TUBElog,"TUBEMODE"
195 call readi(controlcard,'IPRINT',iprint,0)
196 C SHIELD keyword sets if the shielding effect of side-chains is used
197 C 0 denots no shielding is used all peptide are equally despite the
198 C solvent accesible area
199 C 1 the newly introduced function
200 C 2 reseved for further possible developement
201 call readi(controlcard,'SHIELD',shield_mode,0)
202 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
203 write(iout,*) "shield_mode",shield_mode
205 call readi(controlcard,'TORMODE',tor_mode,0)
206 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
207 write(iout,*) "torsional and valence angle mode",tor_mode
208 call readi(controlcard,'MAXGEN',maxgen,10000)
209 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
210 call readi(controlcard,"KDIAG",kdiag,0)
211 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
212 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
213 & write (iout,*) "RESCALE_MODE",rescale_mode
214 split_ene=index(controlcard,'SPLIT_ENE').gt.0
215 if (index(controlcard,'REGULAR').gt.0.0D0) then
216 call reada(controlcard,'WEIDIS',weidis,0.1D0)
221 if (index(controlcard,"CASC").gt.0) then
223 c else if (index(controlcard,"CAONLY").gt.0) then
225 else if (index(controlcard,"SCONLY").gt.0) then
229 c write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
230 c & " specifying REFSTR, PDBREF or REGULAR."
234 if (index(controlcard,'CHECKGRAD').gt.0) then
236 if (index(controlcard,' CART').gt.0) then
238 elseif (index(controlcard,'CARINT').gt.0) then
243 call reada(controlcard,'DELTA',aincr,1.0d-7)
244 c write (iout,*) "icheckgrad",icheckgrad
245 elseif (index(controlcard,'THREAD').gt.0) then
247 call readi(controlcard,'THREAD',nthread,0)
248 if (nthread.gt.0) then
249 call reada(controlcard,'WEIDIS',weidis,0.1D0)
252 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
253 stop 'Error termination in Read_Control.'
255 else if (index(controlcard,'MCMA').gt.0) then
257 else if (index(controlcard,'MCEE').gt.0) then
259 else if (index(controlcard,'MULTCONF').gt.0) then
261 else if (index(controlcard,'MAP').gt.0) then
263 call readi(controlcard,'MAP',nmap,0)
264 else if (index(controlcard,'CSA').gt.0) then
266 crc else if (index(controlcard,'ZSCORE').gt.0) then
268 crc ZSCORE is rm from UNRES, modecalc=9 is available
271 cfcm else if (index(controlcard,'MCMF').gt.0) then
273 else if (index(controlcard,'SOFTREG').gt.0) then
275 else if (index(controlcard,'CHECK_BOND').gt.0) then
277 else if (index(controlcard,'TEST').gt.0) then
279 else if (index(controlcard,'MD').gt.0) then
281 else if (index(controlcard,'RE ').gt.0) then
285 lmuca=index(controlcard,'MUCA').gt.0
286 call readi(controlcard,'MUCADYN',mucadyn,0)
287 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
288 if (lmuca .and. (me.eq.king .or. .not.out1file ))
290 write (iout,*) 'MUCADYN=',mucadyn
291 write (iout,*) 'MUCASMOOTH=',muca_smooth
294 iscode=index(controlcard,'ONE_LETTER')
295 indphi=index(controlcard,'PHI')
296 indback=index(controlcard,'BACK')
297 iranconf=index(controlcard,'RAND_CONF')
298 i2ndstr=index(controlcard,'USE_SEC_PRED')
299 gradout=index(controlcard,'GRADOUT').gt.0
300 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
301 C DISTCHAINMAX become obsolete for periodic boundry condition
302 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
303 C Reading the dimensions of box in x,y,z coordinates
304 call reada(controlcard,'BOXX',boxxsize,100.0d0)
305 call reada(controlcard,'BOXY',boxysize,100.0d0)
306 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
307 write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
308 c Cutoff range for interactions
309 call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
310 call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
311 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
312 write (iout,*) "Cutoff on interactions",r_cut_int
314 & "Cutoff in switching short and long range interactions in RESPA",
316 write (iout,*) "lambda in switch function",rlamb
317 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
318 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
319 if (lipthick.gt.0.0d0) then
320 bordliptop=(boxzsize+lipthick)/2.0
321 bordlipbot=bordliptop-lipthick
323 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
324 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
325 buflipbot=bordlipbot+lipbufthick
326 bufliptop=bordliptop-lipbufthick
327 if ((lipbufthick*2.0d0).gt.lipthick)
328 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
330 write(iout,*) "bordliptop=",bordliptop
331 write(iout,*) "bordlipbot=",bordlipbot
332 write(iout,*) "bufliptop=",bufliptop
333 write(iout,*) "buflipbot=",buflipbot
334 write (iout,*) "SHIELD MODE",shield_mode
335 if (TUBElog.gt.0) then
336 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
337 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
338 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
339 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
340 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
341 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
342 buftubebot=bordtubebot+tubebufthick
343 buftubetop=bordtubetop-tubebufthick
345 if (me.eq.king .or. .not.out1file )
346 & write (iout,*) "DISTCHAINMAX",distchainmax
348 if(me.eq.king.or..not.out1file)
349 & write (iout,'(2a)') diagmeth(kdiag),
350 & ' routine used to diagonalize matrices.'
353 c--------------------------------------------------------------------------
354 subroutine read_REMDpar
360 include 'COMMON.IOUNITS'
361 include 'COMMON.TIME1'
364 include 'COMMON.LANGEVIN'
367 include 'COMMON.LANGEVIN.lang0.5diag'
369 include 'COMMON.LANGEVIN.lang0'
372 include 'COMMON.INTERACT'
373 include 'COMMON.NAMES'
375 include 'COMMON.REMD'
376 include 'COMMON.CONTROL'
377 include 'COMMON.SETUP'
379 character*320 controlcard
380 character*3200 controlcard1
381 integer iremd_m_total,i
383 if(me.eq.king.or..not.out1file)
384 & write (iout,*) "REMD setup"
386 call card_concat(controlcard)
387 call readi(controlcard,"NREP",nrep,3)
388 call readi(controlcard,"NSTEX",nstex,1000)
389 call reada(controlcard,"RETMIN",retmin,10.0d0)
390 call reada(controlcard,"RETMAX",retmax,1000.0d0)
391 mremdsync=(index(controlcard,'SYNC').gt.0)
392 call readi(controlcard,"NSYN",i_sync_step,100)
393 restart1file=(index(controlcard,'REST1FILE').gt.0)
394 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
395 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
396 if(max_cache_traj_use.gt.max_cache_traj)
397 & max_cache_traj_use=max_cache_traj
398 if(me.eq.king.or..not.out1file) then
399 cd if (traj1file) then
400 crc caching is in testing - NTWX is not ignored
401 cd write (iout,*) "NTWX value is ignored"
402 cd write (iout,*) " trajectory is stored to one file by master"
403 cd write (iout,*) " before exchange at NSTEX intervals"
405 write (iout,*) "NREP= ",nrep
406 write (iout,*) "NSTEX= ",nstex
407 write (iout,*) "SYNC= ",mremdsync
408 write (iout,*) "NSYN= ",i_sync_step
409 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
412 if (index(controlcard,'TLIST').gt.0) then
414 call card_concat(controlcard1)
415 read(controlcard1,*) (remd_t(i),i=1,nrep)
416 if(me.eq.king.or..not.out1file)
417 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
420 if (index(controlcard,'MLIST').gt.0) then
422 call card_concat(controlcard1)
423 read(controlcard1,*) (remd_m(i),i=1,nrep)
424 if(me.eq.king.or..not.out1file) then
425 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
428 iremd_m_total=iremd_m_total+remd_m(i)
430 write (iout,*) 'Total number of replicas ',iremd_m_total
435 & "Adaptive (PMF-biased) umbrella sampling will be run"
438 if(me.eq.king.or..not.out1file)
439 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
442 c--------------------------------------------------------------------------
443 subroutine read_MDpar
449 include 'COMMON.IOUNITS'
450 include 'COMMON.TIME1'
452 include 'COMMON.QRESTR'
454 include 'COMMON.LANGEVIN'
457 include 'COMMON.LANGEVIN.lang0.5diag'
459 include 'COMMON.LANGEVIN.lang0'
462 include 'COMMON.INTERACT'
463 include 'COMMON.NAMES'
465 include 'COMMON.SETUP'
466 include 'COMMON.CONTROL'
467 include 'COMMON.SPLITELE'
468 include 'COMMON.FFIELD'
470 character*320 controlcard
474 call card_concat(controlcard)
475 call readi(controlcard,"NSTEP",n_timestep,1000000)
476 call readi(controlcard,"NTWE",ntwe,100)
477 call readi(controlcard,"NTWX",ntwx,1000)
478 call reada(controlcard,"DT",d_time,1.0d-1)
479 call reada(controlcard,"DVMAX",dvmax,2.0d1)
480 call reada(controlcard,"DAMAX",damax,1.0d1)
481 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
482 call readi(controlcard,"LANG",lang,0)
483 RESPA = index(controlcard,"RESPA") .gt. 0
484 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
485 ntime_split0=ntime_split
486 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
487 ntime_split0=ntime_split
488 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
489 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
490 rest = index(controlcard,"REST").gt.0
491 tbf = index(controlcard,"TBF").gt.0
492 usampl = index(controlcard,"USAMPL").gt.0
493 scale_umb = index(controlcard,"SCALE_UMB").gt.0
494 adaptive = index(controlcard,"ADAPTIVE").gt.0
495 mdpdb = index(controlcard,"MDPDB").gt.0
496 call reada(controlcard,"T_BATH",t_bath,300.0d0)
497 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
498 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
499 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
500 if (count_reset_moment.eq.0) count_reset_moment=1000000000
501 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
502 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
503 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
504 if (count_reset_vel.eq.0) count_reset_vel=1000000000
505 large = index(controlcard,"LARGE").gt.0
506 print_compon = index(controlcard,"PRINT_COMPON").gt.0
507 rattle = index(controlcard,"RATTLE").gt.0
508 preminim = index(controlcard,"PREMINIM").gt.0
509 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
513 write (iout,*) "PREMINIM ",preminim
515 if (index(controlcard,'CART').gt.0) dccart=.true.
516 write (iout,*) "dccart ",dccart
517 write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
518 if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
520 c if performing umbrella sampling, fragments constrained are read from the fragment file
523 write (iout,*) "Umbrella sampling will be run"
524 if (scale_umb.and.adaptive) then
525 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
526 write (iout,*) "Select one of those and re-run the job."
529 if (scale_umb) write (iout,*)
530 &"Umbrella-restraint force constants will be scaled by temperature"
534 if(me.eq.king.or..not.out1file) then
536 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
538 write (iout,'(a)') "The units are:"
539 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
540 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
541 & " acceleration: angstrom/(48.9 fs)**2"
542 write (iout,'(a)') "energy: kcal/mol, temperature: K"
544 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
545 write (iout,'(a60,f10.5,a)')
546 & "Initial time step of numerical integration:",d_time,
548 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
550 write (iout,'(2a,i4,a)')
551 & "A-MTS algorithm used; initial time step for fast-varying",
552 & " short-range forces split into",ntime_split," steps."
553 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
554 & r_cut_respa," lambda",rlamb
556 write (iout,'(2a,f10.5)')
557 & "Maximum acceleration threshold to reduce the time step",
558 & "/increase split number:",damax
559 write (iout,'(2a,f10.5)')
560 & "Maximum predicted energy drift to reduce the timestep",
561 & "/increase split number:",edriftmax
562 write (iout,'(a60,f10.5)')
563 & "Maximum velocity threshold to reduce velocities:",dvmax
564 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
565 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
566 if (rattle) write (iout,'(a60)')
567 & "Rattle algorithm used to constrain the virtual bonds"
571 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
572 call reada(controlcard,"RWAT",rwat,1.4d0)
573 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
574 surfarea=index(controlcard,"SURFAREA").gt.0
575 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
576 if(me.eq.king.or..not.out1file)then
577 write (iout,'(/a,$)') "Langevin dynamics calculation"
580 & " with direct integration of Langevin equations"
581 else if (lang.eq.2) then
582 write (iout,'(a/)') " with TINKER stochasic MD integrator"
583 else if (lang.eq.3) then
584 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
585 else if (lang.eq.4) then
586 write (iout,'(a/)') " in overdamped mode"
588 write (iout,'(//a,i5)')
589 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
592 write (iout,'(a60,f10.5)') "Temperature:",t_bath
593 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
594 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
595 write (iout,'(a60,f10.5)')
596 & "Scaling factor of the friction forces:",scal_fric
597 if (surfarea) write (iout,'(2a,i10,a)')
598 & "Friction coefficients will be scaled by solvent-accessible",
599 & " surface area every",reset_fricmat," steps."
601 c Calculate friction coefficients and bounds of stochastic forces
602 eta=6*pi*cPoise*etawat
603 if(me.eq.king.or..not.out1file)
604 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
606 gamp=scal_fric*(pstok+rwat)*eta
607 stdfp=dsqrt(2*Rb*t_bath/d_time)
609 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
610 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
612 if(me.eq.king.or..not.out1file)then
613 write (iout,'(/2a/)')
614 & "Radii of site types and friction coefficients and std's of",
615 & " stochastic forces of fully exposed sites"
616 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
618 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
619 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
623 if(me.eq.king.or..not.out1file)then
624 write (iout,'(a)') "Berendsen bath calculation"
625 write (iout,'(a60,f10.5)') "Temperature:",t_bath
626 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
628 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
629 & count_reset_moment," steps"
631 & write (iout,'(a,i10,a)')
632 & "Velocities will be reset at random every",count_reset_vel,
636 if(me.eq.king.or..not.out1file)
637 & write (iout,'(a31)') "Microcanonical mode calculation"
639 if(me.eq.king.or..not.out1file)then
640 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
642 write(iout,*) "MD running with constraints."
643 write(iout,*) "Equilibration time ", eq_time, " mtus."
644 write(iout,*) "Constraining ", nfrag," fragments."
645 write(iout,*) "Length of each fragment, weight and q0:"
647 write (iout,*) "Set of restraints #",iset
649 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
650 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
652 write(iout,*) "constraints between ", npair, "fragments."
653 write(iout,*) "constraint pairs, weights and q0:"
655 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
656 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
658 write(iout,*) "angle constraints within ", nfrag_back,
659 & "backbone fragments."
661 write(iout,*) "fragment, weights, q0:"
663 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
664 & ifrag_back(2,i,iset),
665 & wfrag_back(1,i,iset),qin_back(1,i,iset),
666 & wfrag_back(2,i,iset),qin_back(2,i,iset),
667 & wfrag_back(3,i,iset),qin_back(3,i,iset)
670 write(iout,*) "fragment, weights:"
672 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
673 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
674 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
678 iset=mod(kolor,nset)+1
681 if(me.eq.king.or..not.out1file)
682 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
685 c------------------------------------------------------------------------------
688 C Read molecular data.
694 integer error_msg,ierror,ierr,ierrcode
696 include 'COMMON.IOUNITS'
699 include 'COMMON.INTERACT'
700 include 'COMMON.LOCAL'
701 include 'COMMON.NAMES'
702 include 'COMMON.CHAIN'
703 include 'COMMON.FFIELD'
704 include 'COMMON.SBRIDGE'
705 include 'COMMON.HEADER'
706 include 'COMMON.CONTROL'
707 include 'COMMON.SAXS'
708 include 'COMMON.DBASE'
709 include 'COMMON.THREAD'
710 include 'COMMON.CONTACTS'
711 include 'COMMON.TORCNSTR'
712 include 'COMMON.TIME1'
713 include 'COMMON.BOUNDS'
715 include 'COMMON.SETUP'
716 include 'COMMON.SHIELD'
717 character*4 sequence(maxres)
719 double precision x(maxvar)
720 character*256 pdbfile
721 character*400 weightcard
722 character*80 weightcard_t,ucase
723 integer itype_pdb(maxres)
724 common /pizda/ itype_pdb
725 logical seq_comp,fail
726 double precision energia(0:n_ene)
727 double precision secprob(3,maxdih_constr)
729 double precision phihel,phibet,sigmahel,sigmabet
730 integer iti,nsi,maxsi
734 integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2
735 double precision sumv
737 C Read PDB structure if applicable
739 if (indpdb.gt.0 .or. pdbref) then
740 read(inp,'(a)') pdbfile
741 if(me.eq.king.or..not.out1file)
742 & write (iout,'(2a)') 'PDB data will be read from file ',
743 & pdbfile(:ilen(pdbfile))
744 open(ipdbin,file=pdbfile,status='old',err=33)
746 33 write (iout,'(a)') 'Error opening PDB file.'
757 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
758 & (crefjlee(j,i+nres),j=1,3)
761 if(me.eq.king.or..not.out1file)
762 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
763 & ' nstart_sup=',nstart_sup
765 itype_pdb(i)=itype(i)
769 nct=nstart_sup+nsup-1
770 call contact(.false.,ncont_ref,icont_ref,co)
773 if(me.eq.king.or..not.out1file)
774 & write(iout,*)'Adding sidechains'
778 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
781 do while (fail.and.nsi.le.maxsi)
782 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
785 if(fail) write(iout,*)'Adding sidechain failed for res ',
786 & i,' after ',nsi,' trials'
791 if (indpdb.eq.0) then
792 C Read sequence if not taken from the pdb file.
794 c print *,'nres=',nres
795 if (iscode.gt.0) then
796 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
798 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
800 C Convert sequence to numeric code
802 itype(i)=rescode(i,sequence(i),iscode)
806 c print '(20i4)',(itype(i),i=1,nres)
809 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
811 if (itype(i).eq.ntyp1) then
815 else if (iabs(itype(i+1)).ne.20) then
817 else if (iabs(itype(i)).ne.20) then
824 if(me.eq.king.or..not.out1file)then
825 write (iout,*) "ITEL"
827 write (iout,*) i,itype(i),itel(i)
829 print *,'Call Read_Bridge.'
833 cd print *,'NNT=',NNT,' NCT=',NCT
834 call seq2chains(nres,itype,nchain,chain_length,chain_border,
837 chain_border1(2,1)=chain_border(2,1)+1
839 chain_border1(1,i)=chain_border(1,i)-1
840 chain_border1(2,i)=chain_border(2,i)+1
842 chain_border1(1,nchain)=chain_border(1,nchain)-1
843 chain_border1(2,nchain)=nres
844 write(iout,*) "nres",nres," nchain",nchain
846 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
847 & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
849 call chain_symmetry(nchain,nres,itype,chain_border,
850 & chain_length,npermchain,tabpermchain)
852 c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
855 write(iout,*) "residue permutations"
857 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
860 if (itype(1).eq.ntyp1) nnt=2
861 if (itype(nres).eq.ntyp1) nct=nct-1
862 write (iout,*) "nnt",nnt," nct",nct
865 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
866 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
868 print*, 'init_dfa_vars finished!'
870 print*, 'read_dfa_info finished!'
874 if(me.eq.king.or..not.out1file)
875 & write (iout,'(a,i3)') 'nsup=',nsup
877 if (nsup.le.(nct-nnt+1)) then
878 do i=0,nct-nnt+1-nsup
879 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
885 & 'Error - sequences to be superposed do not match.'
888 do i=0,nsup-(nct-nnt+1)
889 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
891 nstart_sup=nstart_sup+i
897 & 'Error - sequences to be superposed do not match.'
900 if (nsup.eq.0) nsup=nct-nnt
901 if (nstart_sup.eq.0) nstart_sup=nnt
902 if (nstart_seq.eq.0) nstart_seq=nnt
903 if(me.eq.king.or..not.out1file)
904 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
905 & ' nstart_seq=',nstart_seq
908 C 8/13/98 Set limits to generating the dihedral angles
913 read (inp,*) ndih_constr
914 write (iout,*) "ndih_constr",ndih_constr
915 if (ndih_constr.gt.0) then
918 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
921 if(me.eq.king.or..not.out1file)then
924 & 'There are',ndih_constr,' restraints on gamma angles.'
926 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
933 phi0(i)=deg2rad*phi0(i)
934 drange(i)=deg2rad*drange(i)
936 C if(me.eq.king.or..not.out1file)
937 C & write (iout,*) 'FTORS',ftors
940 phibound(1,ii) = phi0(i)-drange(i)
941 phibound(2,ii) = phi0(i)+drange(i)
944 if (me.eq.king .or. .not.out1file) then
946 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
948 write (iout,'(a3,i5,2f10.1)')
949 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
954 else if (ndih_constr.lt.0) then
956 call card_concat(weightcard)
957 call reada(weightcard,"PHIHEL",phihel,50.0D0)
958 call reada(weightcard,"PHIBET",phibet,180.0D0)
959 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
960 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
961 call reada(weightcard,"WDIHC",wdihc,0.591D0)
962 write (iout,*) "Weight of dihedral angle restraints",wdihc
963 read(inp,'(9x,3f7.3)')
964 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
965 write (iout,*) "The secprob array"
967 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
971 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
972 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
973 ndih_constr=ndih_constr+1
974 idih_constr(ndih_constr)=i
977 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
978 sumv=sumv+vpsipred(j,ndih_constr)
981 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
983 phibound(1,ndih_constr)=phihel*deg2rad
984 phibound(2,ndih_constr)=phibet*deg2rad
985 sdihed(1,ndih_constr)=sigmahel*deg2rad
986 sdihed(2,ndih_constr)=sigmabet*deg2rad
990 if(me.eq.king.or..not.out1file)then
993 & 'There are',ndih_constr,
994 & ' bimodal restraints on gamma angles.'
996 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
997 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
998 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
999 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1000 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1001 & (vpsipred(j,i),j=1,3)
1008 C first setting the theta boundaries to 0 to pi
1009 C this mean that there is no energy penalty for any angle occuring this can be applied
1010 C for generate random conformation but is not implemented in this way
1013 C thetabound(2,i)=pi
1015 C begin reading theta constrains this is quartic constrains allowing to
1016 C have smooth second derivative
1017 if (with_theta_constr) then
1018 C with_theta_constr is keyword allowing for occurance of theta constrains
1019 read (inp,*) ntheta_constr
1020 C ntheta_constr is the number of theta constrains
1021 if (ntheta_constr.gt.0) then
1022 C read (inp,*) ftors
1023 read (inp,*) (itheta_constr(i),theta_constr0(i),
1024 & theta_drange(i),for_thet_constr(i),
1025 & i=1,ntheta_constr)
1026 C the above code reads from 1 to ntheta_constr
1027 C itheta_constr(i) residue i for which is theta_constr
1028 C theta_constr0 the global minimum value
1029 C theta_drange is range for which there is no energy penalty
1030 C for_thet_constr is the force constant for quartic energy penalty
1032 if(me.eq.king.or..not.out1file)then
1034 & 'There are',ntheta_constr,' constraints on phi angles.'
1035 do i=1,ntheta_constr
1036 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1038 & for_thet_constr(i)
1041 do i=1,ntheta_constr
1042 theta_constr0(i)=deg2rad*theta_constr0(i)
1043 theta_drange(i)=deg2rad*theta_drange(i)
1045 C if(me.eq.king.or..not.out1file)
1046 C & write (iout,*) 'FTORS',ftors
1047 C do i=1,ntheta_constr
1048 C ii = itheta_constr(i)
1049 C thetabound(1,ii) = phi0(i)-drange(i)
1050 C thetabound(2,ii) = phi0(i)+drange(i)
1052 endif ! ntheta_constr.gt.0
1053 endif! with_theta_constr
1055 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1056 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1057 c--- Zscore rms -------
1058 if (nz_start.eq.0) nz_start=nnt
1059 if (nz_end.eq.0 .and. nsup.gt.0) then
1061 else if (nz_end.eq.0) then
1064 if(me.eq.king.or..not.out1file)then
1065 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1066 write (iout,*) 'IZ_SC=',iz_sc
1068 c----------------------
1071 if (.not.pdbref) then
1072 call read_angles(inp,*38)
1075 38 write (iout,'(a)') 'Error reading reference structure.'
1077 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1078 stop 'Error reading reference structure'
1080 39 call chainbuild_extconf
1082 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1091 call contact(.true.,ncont_ref,icont_ref,co)
1095 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1097 if (constr_dist.gt.0) call read_dist_constr
1098 write (iout,*) "After read_dist_constr nhpb",nhpb
1099 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1101 call NMRpeak_partition
1102 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1104 if(me.eq.king.or..not.out1file)
1105 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1108 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1110 if(me.eq.king.or..not.out1file)
1111 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1112 & icont_ref(1,i),' ',
1113 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1116 write (iout,*) "calling read_saxs_consrtr",nsaxs
1117 if (nsaxs.gt.0) call read_saxs_constr
1119 if (constr_homology.gt.0) then
1120 call read_constr_homology
1121 if (indpdb.gt.0 .or. pdbref) then
1124 c(j,i)=crefjlee(j,i)
1125 cref(j,i)=crefjlee(j,i)
1130 write (iout,*) "sc_loc_geom: Array C"
1132 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1133 & (c(j,i+nres),j=1,3)
1135 write (iout,*) "Array Cref"
1137 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1138 & (cref(j,i+nres),j=1,3)
1141 call int_from_cart1(.false.)
1142 call sc_loc_geom(.false.)
1144 thetaref(i)=theta(i)
1149 dc(j,i)=c(j,i+1)-c(j,i)
1150 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1155 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1156 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1165 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1166 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1167 & modecalc.ne.10) then
1168 C If input structure hasn't been supplied from the PDB file read or generate
1170 if (iranconf.eq.0 .and. .not. extconf) then
1171 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1172 & write (iout,'(a)') 'Initial geometry will be read in.'
1174 read(inp,'(8f10.5)',end=36,err=36)
1175 & ((c(l,k),l=1,3),k=1,nres),
1176 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1177 write (iout,*) "Exit READ_CART"
1178 c write (iout,'(8f10.5)')
1179 c & ((c(l,k),l=1,3),k=1,nres),
1180 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1184 dc(j,i)=c(j,i+1)-c(j,i)
1185 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1189 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1191 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1192 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1197 c dc_norm(j,i+nres)=0.0d0
1201 call int_from_cart1(.true.)
1202 write (iout,*) "Finish INT_TO_CART"
1203 c write (iout,*) "DC"
1205 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1206 c & (dc(j,i+nres),j=1,3)
1212 call read_angles(inp,*36)
1214 call chainbuild_extconf
1217 36 write (iout,'(a)') 'Error reading angle file.'
1219 call mpi_finalize( MPI_COMM_WORLD,IERR )
1221 stop 'Error reading angle file.'
1223 else if (extconf) then
1224 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1225 & write (iout,'(a)') 'Extended chain initial geometry.'
1227 theta(i)=90d0*deg2rad
1230 phi(i)=180d0*deg2rad
1233 alph(i)=110d0*deg2rad
1236 omeg(i)=-120d0*deg2rad
1237 if (itype(i).le.0) omeg(i)=-omeg(i)
1240 call chainbuild_extconf
1242 if(me.eq.king.or..not.out1file)
1243 & write (iout,'(a)') 'Random-generated initial geometry.'
1246 if (me.eq.king .or. fg_rank.eq.0 .and. (
1247 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1251 call gen_rand_conf(itmp,*30)
1253 30 write (iout,*) 'Failed to generate random conformation',
1254 & ', itrial=',itrial
1255 write (*,*) 'Processor:',me,
1256 & ' Failed to generate random conformation',
1266 write (iout,'(a,i3,a)') 'Processor:',me,
1267 & ' error in generating random conformation.'
1268 write (*,'(a,i3,a)') 'Processor:',me,
1269 & ' error in generating random conformation.'
1272 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1277 & ' error in generating random conformation.'
1282 elseif (modecalc.eq.4) then
1283 read (inp,'(a)') intinname
1284 open (intin,file=intinname,status='old',err=333)
1285 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1286 & write (iout,'(a)') 'intinname',intinname
1287 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1289 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1291 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1293 stop 'Error opening angle file.'
1297 C Generate distance constraints, if the PDB structure is to be regularized.
1298 if (nthread.gt.0) then
1299 call read_threadbase
1302 if (me.eq.king .or. .not. out1file)
1304 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1305 write (iout,'(/a,i3,a)')
1306 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1307 write (iout,'(20i4)') (iss(i),i=1,ns)
1309 write(iout,*)"Running with dynamic disulfide-bond formation"
1311 write (iout,'(/a/)') 'Pre-formed links are:'
1317 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1318 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1324 if (ns.gt.0.and.dyn_ss) then
1328 forcon(i-nss)=forcon(i)
1335 dyn_ss_mask(iss(i))=.true.
1340 c write (iout,*) "DC"
1342 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1343 c & (dc(j,i+nres),j=1,3)
1345 if (i2ndstr.gt.0) call secstrp2dihc
1346 c call geom_to_var(nvar,x)
1347 c call etotal(energia(0))
1348 c call enerprint(energia(0))
1349 c call briefout(0,etot)
1351 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1352 cd write (iout,'(a)') 'Variable list:'
1353 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1355 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1356 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1357 & 'Processor',myrank,': end reading molecular data.'
1362 c--------------------------------------------------------------------------
1363 logical function seq_comp(itypea,itypeb,length)
1365 integer length,itypea(length),itypeb(length)
1368 if (itypea(i).ne.itypeb(i)) then
1376 c-----------------------------------------------------------------------------
1377 subroutine read_bridge
1378 C Read information about disulfide bridges.
1380 include 'DIMENSIONS'
1385 include 'COMMON.IOUNITS'
1386 include 'COMMON.GEO'
1387 include 'COMMON.VAR'
1388 include 'COMMON.INTERACT'
1389 include 'COMMON.LOCAL'
1390 include 'COMMON.NAMES'
1391 include 'COMMON.CHAIN'
1392 include 'COMMON.FFIELD'
1393 include 'COMMON.SBRIDGE'
1394 include 'COMMON.HEADER'
1395 include 'COMMON.CONTROL'
1396 include 'COMMON.DBASE'
1397 include 'COMMON.THREAD'
1398 include 'COMMON.TIME1'
1399 include 'COMMON.SETUP'
1401 C Read bridging residues.
1402 read (inp,*) ns,(iss(i),i=1,ns)
1404 if(me.eq.king.or..not.out1file)
1405 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1406 C Check whether the specified bridging residues are cystines.
1408 if (itype(iss(i)).ne.1) then
1409 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1410 & 'Do you REALLY think that the residue ',
1411 & restyp(itype(iss(i))),i,
1412 & ' can form a disulfide bridge?!!!'
1413 write (*,'(2a,i3,a)')
1414 & 'Do you REALLY think that the residue ',
1415 & restyp(itype(iss(i))),i,
1416 & ' can form a disulfide bridge?!!!'
1418 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1423 C Read preformed bridges.
1425 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1427 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1430 C Check if the residues involved in bridges are in the specified list of
1431 C bridging residues.
1434 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1435 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1436 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1437 & ' contains residues present in other pairs.'
1438 write (*,'(a,i3,a)') 'Disulfide pair',i,
1439 & ' contains residues present in other pairs.'
1441 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1447 if (ihpb(i).eq.iss(j)) goto 10
1449 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1452 if (jhpb(i).eq.iss(j)) goto 20
1454 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1460 ihpb(i)=ihpb(i)+nres
1461 jhpb(i)=jhpb(i)+nres
1467 c----------------------------------------------------------------------------
1468 subroutine read_x(kanal,*)
1470 include 'DIMENSIONS'
1471 include 'COMMON.GEO'
1472 include 'COMMON.VAR'
1473 include 'COMMON.CHAIN'
1474 include 'COMMON.IOUNITS'
1475 include 'COMMON.CONTROL'
1476 include 'COMMON.LOCAL'
1477 include 'COMMON.INTERACT'
1478 integer i,j,k,l,kanal
1479 c Read coordinates from input
1481 read(kanal,'(8f10.5)',end=10,err=10)
1482 & ((c(l,k),l=1,3),k=1,nres),
1483 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1486 c(j,2*nres)=c(j,nres)
1488 call int_from_cart1(.false.)
1491 dc(j,i)=c(j,i+1)-c(j,i)
1492 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1496 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1498 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1499 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1507 c----------------------------------------------------------------------------
1508 subroutine read_threadbase
1510 include 'DIMENSIONS'
1511 include 'COMMON.IOUNITS'
1512 include 'COMMON.GEO'
1513 include 'COMMON.VAR'
1514 include 'COMMON.INTERACT'
1515 include 'COMMON.LOCAL'
1516 include 'COMMON.NAMES'
1517 include 'COMMON.CHAIN'
1518 include 'COMMON.FFIELD'
1519 include 'COMMON.SBRIDGE'
1520 include 'COMMON.HEADER'
1521 include 'COMMON.CONTROL'
1522 include 'COMMON.DBASE'
1523 include 'COMMON.THREAD'
1524 include 'COMMON.TIME1'
1526 double precision dist
1527 C Read pattern database for threading.
1528 read (icbase,*) nseq
1530 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1531 & nres_base(2,i),nres_base(3,i)
1532 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1534 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1535 c & nres_base(2,i),nres_base(3,i)
1536 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1540 if (weidis.eq.0.0D0) weidis=0.1D0
1549 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1550 write (iout,'(a,i5)') 'nexcl: ',nexcl
1551 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1554 c------------------------------------------------------------------------------
1555 subroutine setup_var
1557 include 'DIMENSIONS'
1558 include 'COMMON.IOUNITS'
1559 include 'COMMON.GEO'
1560 include 'COMMON.VAR'
1561 include 'COMMON.INTERACT'
1562 include 'COMMON.LOCAL'
1563 include 'COMMON.NAMES'
1564 include 'COMMON.CHAIN'
1565 include 'COMMON.FFIELD'
1566 include 'COMMON.SBRIDGE'
1567 include 'COMMON.HEADER'
1568 include 'COMMON.CONTROL'
1569 include 'COMMON.DBASE'
1570 include 'COMMON.THREAD'
1571 include 'COMMON.TIME1'
1573 C Set up variable list.
1579 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1581 ialph(i,1)=nvar+nside
1585 if (indphi.gt.0) then
1587 else if (indback.gt.0) then
1594 c----------------------------------------------------------------------------
1595 subroutine gen_dist_constr
1596 C Generate CA distance constraints.
1598 include 'DIMENSIONS'
1599 include 'COMMON.IOUNITS'
1600 include 'COMMON.GEO'
1601 include 'COMMON.VAR'
1602 include 'COMMON.INTERACT'
1603 include 'COMMON.LOCAL'
1604 include 'COMMON.NAMES'
1605 include 'COMMON.CHAIN'
1606 include 'COMMON.FFIELD'
1607 include 'COMMON.SBRIDGE'
1608 include 'COMMON.HEADER'
1609 include 'COMMON.CONTROL'
1610 include 'COMMON.DBASE'
1611 include 'COMMON.THREAD'
1612 include 'COMMON.TIME1'
1613 integer i,j,itype_pdb(maxres)
1614 common /pizda/ itype_pdb
1615 double precision dist
1617 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1618 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1619 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1621 do i=nstart_sup,nstart_sup+nsup-1
1622 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1623 cd & ' seq_pdb', restyp(itype_pdb(i))
1624 do j=i+2,nstart_sup+nsup-1
1626 ihpb(nhpb)=i+nstart_seq-nstart_sup
1627 jhpb(nhpb)=j+nstart_seq-nstart_sup
1629 dhpb(nhpb)=dist(i,j)
1632 cd write (iout,'(a)') 'Distance constraints:'
1637 cd if (ii.gt.nres) then
1642 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1643 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1644 cd & dhpb(i),forcon(i)
1648 c----------------------------------------------------------------------------
1651 include 'DIMENSIONS'
1652 include 'COMMON.MAP'
1653 include 'COMMON.IOUNITS'
1655 character*3 angid(4) /'THE','PHI','ALP','OME'/
1656 character*80 mapcard,ucase
1658 read (inp,'(a)') mapcard
1659 mapcard=ucase(mapcard)
1660 if (index(mapcard,'PHI').gt.0) then
1662 else if (index(mapcard,'THE').gt.0) then
1664 else if (index(mapcard,'ALP').gt.0) then
1666 else if (index(mapcard,'OME').gt.0) then
1669 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1670 stop 'Error - illegal variable spec in MAP card.'
1672 call readi (mapcard,'RES1',res1(imap),0)
1673 call readi (mapcard,'RES2',res2(imap),0)
1674 if (res1(imap).eq.0) then
1675 res1(imap)=res2(imap)
1676 else if (res2(imap).eq.0) then
1677 res2(imap)=res1(imap)
1679 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1681 & 'Error - illegal definition of variable group in MAP.'
1682 stop 'Error - illegal definition of variable group in MAP.'
1684 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1685 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1686 call readi(mapcard,'NSTEP',nstep(imap),0)
1687 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1689 & 'Illegal boundary and/or step size specification in MAP.'
1690 stop 'Illegal boundary and/or step size specification in MAP.'
1695 c----------------------------------------------------------------------------
1698 include 'DIMENSIONS'
1699 include 'COMMON.IOUNITS'
1700 include 'COMMON.GEO'
1701 include 'COMMON.CSA'
1702 include 'COMMON.BANK'
1703 include 'COMMON.CONTROL'
1705 character*620 mcmcard
1706 call card_concat(mcmcard)
1708 call readi(mcmcard,'NCONF',nconf,50)
1709 call readi(mcmcard,'NADD',nadd,0)
1710 call readi(mcmcard,'JSTART',jstart,1)
1711 call readi(mcmcard,'JEND',jend,1)
1712 call readi(mcmcard,'NSTMAX',nstmax,500000)
1713 call readi(mcmcard,'N0',n0,1)
1714 call readi(mcmcard,'N1',n1,6)
1715 call readi(mcmcard,'N2',n2,4)
1716 call readi(mcmcard,'N3',n3,0)
1717 call readi(mcmcard,'N4',n4,0)
1718 call readi(mcmcard,'N5',n5,0)
1719 call readi(mcmcard,'N6',n6,10)
1720 call readi(mcmcard,'N7',n7,0)
1721 call readi(mcmcard,'N8',n8,0)
1722 call readi(mcmcard,'N9',n9,0)
1723 call readi(mcmcard,'N14',n14,0)
1724 call readi(mcmcard,'N15',n15,0)
1725 call readi(mcmcard,'N16',n16,0)
1726 call readi(mcmcard,'N17',n17,0)
1727 call readi(mcmcard,'N18',n18,0)
1729 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1731 call readi(mcmcard,'NDIFF',ndiff,2)
1732 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1733 call readi(mcmcard,'IS1',is1,1)
1734 call readi(mcmcard,'IS2',is2,8)
1735 call readi(mcmcard,'NRAN0',nran0,4)
1736 call readi(mcmcard,'NRAN1',nran1,2)
1737 call readi(mcmcard,'IRR',irr,1)
1738 call readi(mcmcard,'NSEED',nseed,20)
1739 call readi(mcmcard,'NTOTAL',ntotal,10000)
1740 call reada(mcmcard,'CUT1',cut1,2.0d0)
1741 call reada(mcmcard,'CUT2',cut2,5.0d0)
1742 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1743 call readi(mcmcard,'ICMAX',icmax,3)
1744 call readi(mcmcard,'IRESTART',irestart,0)
1745 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1748 call reada(mcmcard,'DELE',dele,20.0d0)
1749 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1750 call readi(mcmcard,'IREF',iref,0)
1751 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1752 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1753 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1754 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1755 write (iout,*) "NCONF_IN",nconf_in
1758 c----------------------------------------------------------------------------
1761 include 'DIMENSIONS'
1762 include 'COMMON.MCM'
1763 include 'COMMON.MCE'
1764 include 'COMMON.IOUNITS'
1766 character*320 mcmcard
1768 call card_concat(mcmcard)
1769 call readi(mcmcard,'MAXACC',maxacc,100)
1770 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1771 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1772 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1773 call readi(mcmcard,'MAXREPM',maxrepm,200)
1774 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1775 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1776 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1777 call reada(mcmcard,'E_UP',e_up,5.0D0)
1778 call reada(mcmcard,'DELTE',delte,0.1D0)
1779 call readi(mcmcard,'NSWEEP',nsweep,5)
1780 call readi(mcmcard,'NSTEPH',nsteph,0)
1781 call readi(mcmcard,'NSTEPC',nstepc,0)
1782 call reada(mcmcard,'TMIN',tmin,298.0D0)
1783 call reada(mcmcard,'TMAX',tmax,298.0D0)
1784 call readi(mcmcard,'NWINDOW',nwindow,0)
1785 call readi(mcmcard,'PRINT_MC',print_mc,0)
1786 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1787 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1788 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1789 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1790 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1791 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1792 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1793 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1794 if (nwindow.gt.0) then
1795 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1797 winlen(i)=winend(i)-winstart(i)+1
1800 if (tmax.lt.tmin) tmax=tmin
1801 if (tmax.eq.tmin) then
1805 if (nstepc.gt.0 .and. nsteph.gt.0) then
1806 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1807 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1809 C Probabilities of different move types
1810 sumpro_type(0)=0.0D0
1811 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1812 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1813 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1814 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1815 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1816 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1817 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1819 print *,'i',i,' sumprotype',sumpro_type(i)
1820 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1821 print *,'i',i,' sumprotype',sumpro_type(i)
1825 c----------------------------------------------------------------------------
1826 subroutine read_minim
1828 include 'DIMENSIONS'
1829 include 'COMMON.MINIM'
1830 include 'COMMON.IOUNITS'
1832 character*320 minimcard
1833 call card_concat(minimcard)
1834 call readi(minimcard,'MAXMIN',maxmin,2000)
1835 call readi(minimcard,'MAXFUN',maxfun,5000)
1836 call readi(minimcard,'MINMIN',minmin,maxmin)
1837 call readi(minimcard,'MINFUN',minfun,maxmin)
1838 call reada(minimcard,'TOLF',tolf,1.0D-2)
1839 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1840 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1841 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1842 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1843 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1844 & 'Options in energy minimization:'
1845 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1846 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1847 & 'MinMin:',MinMin,' MinFun:',MinFun,
1848 & ' TolF:',TolF,' RTolF:',RTolF
1851 c----------------------------------------------------------------------------
1852 subroutine read_angles(kanal,*)
1854 include 'DIMENSIONS'
1855 include 'COMMON.GEO'
1856 include 'COMMON.VAR'
1857 include 'COMMON.CHAIN'
1858 include 'COMMON.IOUNITS'
1859 include 'COMMON.CONTROL'
1861 c Read angles from input
1863 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1864 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1865 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1866 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1869 c 9/7/01 avoid 180 deg valence angle
1870 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1872 theta(i)=deg2rad*theta(i)
1873 phi(i)=deg2rad*phi(i)
1874 alph(i)=deg2rad*alph(i)
1875 omeg(i)=deg2rad*omeg(i)
1880 c----------------------------------------------------------------------------
1881 subroutine reada(rekord,lancuch,wartosc,default)
1883 character*(*) rekord,lancuch
1884 double precision wartosc,default
1887 iread=index(rekord,lancuch)
1888 if (iread.eq.0) then
1892 iread=iread+ilen(lancuch)+1
1893 read (rekord(iread:),*,err=10,end=10) wartosc
1898 c----------------------------------------------------------------------------
1899 subroutine readi(rekord,lancuch,wartosc,default)
1901 character*(*) rekord,lancuch
1902 integer wartosc,default
1905 iread=index(rekord,lancuch)
1906 if (iread.eq.0) then
1910 iread=iread+ilen(lancuch)+1
1911 read (rekord(iread:),*,err=10,end=10) wartosc
1916 c----------------------------------------------------------------------------
1917 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1920 integer tablica(dim),default
1921 character*(*) rekord,lancuch
1928 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1929 if (iread.eq.0) return
1930 iread=iread+ilen(lancuch)+1
1931 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1934 c----------------------------------------------------------------------------
1935 subroutine multreada(rekord,lancuch,tablica,dim,default)
1938 double precision tablica(dim),default
1939 character*(*) rekord,lancuch
1946 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1947 if (iread.eq.0) return
1948 iread=iread+ilen(lancuch)+1
1949 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1952 c----------------------------------------------------------------------------
1953 subroutine openunits
1955 include 'DIMENSIONS'
1959 character*16 form,nodename
1962 include 'COMMON.SETUP'
1963 include 'COMMON.IOUNITS'
1965 include 'COMMON.CONTROL'
1966 integer lenpre,lenpot,ilen,lentmp,npos
1968 character*3 out1file_text,ucase
1973 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1974 call getenv_loc("PREFIX",prefix)
1976 call getenv_loc("POT",pot)
1977 call getenv_loc("DIRTMP",tmpdir)
1978 call getenv_loc("CURDIR",curdir)
1979 call getenv_loc("OUT1FILE",out1file_text)
1980 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1981 out1file_text=ucase(out1file_text)
1982 if (out1file_text(1:1).eq."Y") then
1985 out1file=fg_rank.gt.0
1990 if (lentmp.gt.0) then
1991 write (*,'(80(1h!))')
1992 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1993 write (*,'(80(1h!))')
1994 write (*,*)"All output files will be on node /tmp directory."
1996 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1997 if (me.eq.king) then
1998 write (*,*) "The master node is ",nodename
1999 else if (fg_rank.eq.0) then
2000 write (*,*) "I am the CG slave node ",nodename
2002 write (*,*) "I am the FG slave node ",nodename
2005 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2006 lenpre = lentmp+lenpre+1
2008 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2009 C Get the names and open the input files
2010 #if defined(WINIFL) || defined(WINPGI)
2011 open(1,file=pref_orig(:ilen(pref_orig))//
2012 & '.inp',status='old',readonly,shared)
2013 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2014 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2015 C Get parameter filenames and open the parameter files.
2016 call getenv_loc('BONDPAR',bondname)
2017 open (ibond,file=bondname,status='old',readonly,shared)
2018 call getenv_loc('THETPAR',thetname)
2019 open (ithep,file=thetname,status='old',readonly,shared)
2020 call getenv_loc('ROTPAR',rotname)
2021 open (irotam,file=rotname,status='old',readonly,shared)
2022 call getenv_loc('TORPAR',torname)
2023 open (itorp,file=torname,status='old',readonly,shared)
2024 call getenv_loc('TORDPAR',tordname)
2025 open (itordp,file=tordname,status='old',readonly,shared)
2026 call getenv_loc('FOURIER',fouriername)
2027 open (ifourier,file=fouriername,status='old',readonly,shared)
2028 call getenv_loc('ELEPAR',elename)
2029 open (ielep,file=elename,status='old',readonly,shared)
2030 call getenv_loc('SIDEPAR',sidename)
2031 open (isidep,file=sidename,status='old',readonly,shared)
2032 call getenv_loc('LIPTRANPAR',liptranname)
2033 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2034 #elif (defined CRAY) || (defined AIX)
2035 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2037 c print *,"Processor",myrank," opened file 1"
2038 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2039 c print *,"Processor",myrank," opened file 9"
2040 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2041 C Get parameter filenames and open the parameter files.
2042 call getenv_loc('BONDPAR',bondname)
2043 open (ibond,file=bondname,status='old',action='read')
2044 c print *,"Processor",myrank," opened file IBOND"
2045 call getenv_loc('THETPAR',thetname)
2046 open (ithep,file=thetname,status='old',action='read')
2048 call getenv_loc('THETPARPDB',thetname_pdb)
2049 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2051 c print *,"Processor",myrank," opened file ITHEP"
2052 call getenv_loc('ROTPAR',rotname)
2053 open (irotam,file=rotname,status='old',action='read')
2055 call getenv_loc('ROTPARPDB',rotname_pdb)
2056 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2058 c print *,"Processor",myrank," opened file IROTAM"
2059 call getenv_loc('TORPAR',torname)
2060 open (itorp,file=torname,status='old',action='read')
2061 c print *,"Processor",myrank," opened file ITORP"
2062 call getenv_loc('TORDPAR',tordname)
2063 open (itordp,file=tordname,status='old',action='read')
2064 c print *,"Processor",myrank," opened file ITORDP"
2065 call getenv_loc('SCCORPAR',sccorname)
2066 open (isccor,file=sccorname,status='old',action='read')
2067 c print *,"Processor",myrank," opened file ISCCOR"
2068 call getenv_loc('FOURIER',fouriername)
2069 open (ifourier,file=fouriername,status='old',action='read')
2070 c print *,"Processor",myrank," opened file IFOURIER"
2071 call getenv_loc('ELEPAR',elename)
2072 open (ielep,file=elename,status='old',action='read')
2073 c print *,"Processor",myrank," opened file IELEP"
2074 call getenv_loc('SIDEPAR',sidename)
2075 open (isidep,file=sidename,status='old',action='read')
2076 call getenv_loc('LIPTRANPAR',liptranname)
2077 open (iliptranpar,file=liptranname,status='old',action='read')
2078 c print *,"Processor",myrank," opened file ISIDEP"
2079 c print *,"Processor",myrank," opened parameter files"
2081 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2082 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2083 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2084 C Get parameter filenames and open the parameter files.
2085 call getenv_loc('BONDPAR',bondname)
2086 open (ibond,file=bondname,status='old')
2087 call getenv_loc('THETPAR',thetname)
2088 open (ithep,file=thetname,status='old')
2090 call getenv_loc('THETPARPDB',thetname_pdb)
2091 open (ithep_pdb,file=thetname_pdb,status='old')
2093 call getenv_loc('ROTPAR',rotname)
2094 open (irotam,file=rotname,status='old')
2096 call getenv_loc('ROTPARPDB',rotname_pdb)
2097 open (irotam_pdb,file=rotname_pdb,status='old')
2099 call getenv_loc('TORPAR',torname)
2100 open (itorp,file=torname,status='old')
2101 call getenv_loc('TORDPAR',tordname)
2102 open (itordp,file=tordname,status='old')
2103 call getenv_loc('SCCORPAR',sccorname)
2104 open (isccor,file=sccorname,status='old')
2105 call getenv_loc('FOURIER',fouriername)
2106 open (ifourier,file=fouriername,status='old')
2107 call getenv_loc('ELEPAR',elename)
2108 open (ielep,file=elename,status='old')
2109 call getenv_loc('SIDEPAR',sidename)
2110 open (isidep,file=sidename,status='old')
2111 call getenv_loc('LIPTRANPAR',liptranname)
2112 open (iliptranpar,file=liptranname,status='old')
2114 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2116 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2117 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2118 C Get parameter filenames and open the parameter files.
2119 call getenv_loc('BONDPAR',bondname)
2120 open (ibond,file=bondname,status='old',readonly)
2121 call getenv_loc('THETPAR',thetname)
2122 open (ithep,file=thetname,status='old',readonly)
2123 call getenv_loc('ROTPAR',rotname)
2124 open (irotam,file=rotname,status='old',readonly)
2125 call getenv_loc('TORPAR',torname)
2126 open (itorp,file=torname,status='old',readonly)
2127 call getenv_loc('TORDPAR',tordname)
2128 open (itordp,file=tordname,status='old',readonly)
2129 call getenv_loc('SCCORPAR',sccorname)
2130 open (isccor,file=sccorname,status='old',readonly)
2132 call getenv_loc('THETPARPDB',thetname_pdb)
2133 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2135 call getenv_loc('FOURIER',fouriername)
2136 open (ifourier,file=fouriername,status='old',readonly)
2137 call getenv_loc('ELEPAR',elename)
2138 open (ielep,file=elename,status='old',readonly)
2139 call getenv_loc('SIDEPAR',sidename)
2140 open (isidep,file=sidename,status='old',readonly)
2141 call getenv_loc('LIPTRANPAR',liptranname)
2142 open (iliptranpar,file=liptranname,status='old',action='read')
2144 call getenv_loc('ROTPARPDB',rotname_pdb)
2145 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2149 call getenv_loc('TUBEPAR',tubename)
2150 #if defined(WINIFL) || defined(WINPGI)
2151 open (itube,file=tubename,status='old',readonly,shared)
2152 #elif (defined CRAY) || (defined AIX)
2153 open (itube,file=tubename,status='old',action='read')
2155 open (itube,file=tubename,status='old')
2157 open (itube,file=tubename,status='old',readonly)
2162 C 8/9/01 In the newest version SCp interaction constants are read from a file
2163 C Use -DOLDSCP to use hard-coded constants instead.
2165 call getenv_loc('SCPPAR',scpname)
2166 #if defined(WINIFL) || defined(WINPGI)
2167 open (iscpp,file=scpname,status='old',readonly,shared)
2168 #elif (defined CRAY) || (defined AIX)
2169 open (iscpp,file=scpname,status='old',action='read')
2171 open (iscpp,file=scpname,status='old')
2173 open (iscpp,file=scpname,status='old',readonly)
2176 call getenv_loc('PATTERN',patname)
2177 #if defined(WINIFL) || defined(WINPGI)
2178 open (icbase,file=patname,status='old',readonly,shared)
2179 #elif (defined CRAY) || (defined AIX)
2180 open (icbase,file=patname,status='old',action='read')
2182 open (icbase,file=patname,status='old')
2184 open (icbase,file=patname,status='old',readonly)
2187 C Open output file only for CG processes
2188 c print *,"Processor",myrank," fg_rank",fg_rank
2189 if (fg_rank.eq.0) then
2191 if (nodes.eq.1) then
2194 npos = dlog10(dfloat(nodes-1))+1
2196 if (npos.lt.3) npos=3
2197 write (liczba,'(i1)') npos
2198 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2200 write (liczba,form) me
2201 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2202 & liczba(:ilen(liczba))
2203 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2205 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2207 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2208 & liczba(:ilen(liczba))//'.mol2'
2209 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2210 & liczba(:ilen(liczba))//'.stat'
2212 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2213 & //liczba(:ilen(liczba))//'.stat')
2214 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2217 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2218 & liczba(:ilen(liczba))//'.const'
2223 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2224 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2225 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2226 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2227 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2229 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2231 rest2name=prefix(:ilen(prefix))//'.rst'
2233 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2236 #if defined(AIX) || defined(PGI) || defined(CRAY)
2237 if (me.eq.king .or. .not. out1file) then
2238 open(iout,file=outname,status='unknown')
2240 open(iout,file="/dev/null",status="unknown")
2244 if (fg_rank.gt.0) then
2245 write (liczba,'(i3.3)') myrank/nfgtasks
2246 write (ll,'(bz,i3.3)') fg_rank
2247 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2253 open(igeom,file=intname,status='unknown',position='append')
2254 open(ipdb,file=pdbname,status='unknown')
2255 open(imol2,file=mol2name,status='unknown')
2256 open(istat,file=statname,status='unknown',position='append')
2258 c1out open(iout,file=outname,status='unknown')
2261 if (me.eq.king .or. .not.out1file)
2262 & open(iout,file=outname,status='unknown')
2264 if (fg_rank.gt.0) then
2265 write (liczba,'(i3.3)') myrank/nfgtasks
2266 write (ll,'(bz,i3.3)') fg_rank
2267 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2272 open(igeom,file=intname,status='unknown',access='append')
2273 open(ipdb,file=pdbname,status='unknown')
2274 open(imol2,file=mol2name,status='unknown')
2275 open(istat,file=statname,status='unknown',access='append')
2277 c1out open(iout,file=outname,status='unknown')
2280 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2281 csa_seed=prefix(:lenpre)//'.CSA.seed'
2282 csa_history=prefix(:lenpre)//'.CSA.history'
2283 csa_bank=prefix(:lenpre)//'.CSA.bank'
2284 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2285 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2286 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2287 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2288 csa_int=prefix(:lenpre)//'.int'
2289 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2290 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2291 csa_in=prefix(:lenpre)//'.CSA.in'
2292 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2295 write (iout,'(80(1h-))')
2296 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2297 write (iout,'(80(1h-))')
2298 write (iout,*) "Input file : ",
2299 & pref_orig(:ilen(pref_orig))//'.inp'
2300 write (iout,*) "Output file : ",
2301 & outname(:ilen(outname))
2303 write (iout,*) "Sidechain potential file : ",
2304 & sidename(:ilen(sidename))
2306 write (iout,*) "SCp potential file : ",
2307 & scpname(:ilen(scpname))
2309 write (iout,*) "Electrostatic potential file : ",
2310 & elename(:ilen(elename))
2311 write (iout,*) "Cumulant coefficient file : ",
2312 & fouriername(:ilen(fouriername))
2313 write (iout,*) "Torsional parameter file : ",
2314 & torname(:ilen(torname))
2315 write (iout,*) "Double torsional parameter file : ",
2316 & tordname(:ilen(tordname))
2317 write (iout,*) "SCCOR parameter file : ",
2318 & sccorname(:ilen(sccorname))
2319 write (iout,*) "Bond & inertia constant file : ",
2320 & bondname(:ilen(bondname))
2321 write (iout,*) "Bending-torsion parameter file : ",
2322 & thetname(:ilen(thetname))
2323 write (iout,*) "Rotamer parameter file : ",
2324 & rotname(:ilen(rotname))
2325 write (iout,*) "Thetpdb parameter file : ",
2326 & thetname_pdb(:ilen(thetname_pdb))
2327 write (iout,*) "Threading database : ",
2328 & patname(:ilen(patname))
2330 &write (iout,*)" DIRTMP : ",
2332 write (iout,'(80(1h-))')
2336 c----------------------------------------------------------------------------
2337 subroutine card_concat(card)
2338 implicit real*8 (a-h,o-z)
2339 include 'DIMENSIONS'
2340 include 'COMMON.IOUNITS'
2342 character*80 karta,ucase
2344 read (inp,'(a)') karta
2347 do while (karta(80:80).eq.'&')
2348 card=card(:ilen(card)+1)//karta(:79)
2349 read (inp,'(a)') karta
2352 card=card(:ilen(card)+1)//karta
2355 c------------------------------------------------------------------------------
2358 include 'DIMENSIONS'
2359 include 'COMMON.CHAIN'
2360 include 'COMMON.IOUNITS'
2361 include 'COMMON.CONTROL'
2363 include 'COMMON.QRESTR'
2365 include 'COMMON.LAGRANGE.5diag'
2367 include 'COMMON.LAGRANGE'
2370 open(irest2,file=rest2name,status='unknown')
2371 read(irest2,*) totT,EK,potE,totE,t_bath
2374 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2377 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2380 read (irest2,*) iset
2385 c------------------------------------------------------------------------------
2386 subroutine read_fragments
2388 include 'DIMENSIONS'
2392 include 'COMMON.SETUP'
2393 include 'COMMON.CHAIN'
2394 include 'COMMON.IOUNITS'
2396 include 'COMMON.QRESTR'
2397 include 'COMMON.CONTROL'
2399 read(inp,*) nset,nfrag,npair,nfrag_back
2400 loc_qlike=(nfrag_back.lt.0)
2401 nfrag_back=iabs(nfrag_back)
2402 if(me.eq.king.or..not.out1file)
2403 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2404 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2406 read(inp,*) mset(iset)
2408 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2410 if(me.eq.king.or..not.out1file)
2411 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2412 & ifrag(2,i,iset), qinfrag(i,iset)
2415 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2417 if(me.eq.king.or..not.out1file)
2418 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2419 & ipair(2,i,iset), qinpair(i,iset)
2423 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2424 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2425 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2426 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2427 if(me.eq.king.or..not.out1file)
2428 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2429 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2430 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2431 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2435 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2436 & wfrag_back(3,i,iset),
2437 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2438 if(me.eq.king.or..not.out1file)
2439 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2440 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2446 C---------------------------------------------------------------------------
2447 subroutine read_afminp
2449 include 'DIMENSIONS'
2453 include 'COMMON.SETUP'
2454 include 'COMMON.CONTROL'
2455 include 'COMMON.CHAIN'
2456 include 'COMMON.IOUNITS'
2457 include 'COMMON.SBRIDGE'
2458 character*320 afmcard
2460 c print *, "wchodze"
2461 call card_concat(afmcard)
2462 call readi(afmcard,"BEG",afmbeg,0)
2463 call readi(afmcard,"END",afmend,0)
2464 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2465 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2466 print *,'FORCE=' ,forceAFMconst
2467 CCCC NOW PROPERTIES FOR AFM
2470 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2472 distafminit=dsqrt(distafminit)
2473 c print *,'initdist',distafminit
2476 c-------------------------------------------------------------------------------
2477 subroutine read_saxs_constr
2479 include 'DIMENSIONS'
2483 include 'COMMON.SETUP'
2484 include 'COMMON.CONTROL'
2485 include 'COMMON.SAXS'
2486 include 'COMMON.CHAIN'
2487 include 'COMMON.IOUNITS'
2488 include 'COMMON.SBRIDGE'
2489 double precision cm(3),cnorm
2492 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2494 if (saxs_mode.eq.0) then
2495 c SAXS distance distribution
2497 read(inp,*) distsaxs(i),Psaxs(i)
2501 Cnorm = Cnorm + Psaxs(i)
2503 write (iout,*) "Cnorm",Cnorm
2505 Psaxs(i)=Psaxs(i)/Cnorm
2507 write (iout,*) "Normalized distance distribution from SAXS"
2509 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2513 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2515 write (iout,*) "Wsaxs0",Wsaxs0
2519 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2526 cm(j)=cm(j)+Csaxs(j,i)
2534 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2537 write (iout,*) "SAXS sphere coordinates"
2539 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2545 c-------------------------------------------------------------------------------
2546 subroutine read_dist_constr
2548 include 'DIMENSIONS'
2552 include 'COMMON.SETUP'
2553 include 'COMMON.CONTROL'
2554 include 'COMMON.CHAIN'
2555 include 'COMMON.IOUNITS'
2556 include 'COMMON.SBRIDGE'
2557 include 'COMMON.INTERACT'
2558 integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2559 integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2560 double precision wfrag_(100),wpair_(1000)
2561 double precision ddjk,dist,dist_cut,fordepthmax
2562 character*5000 controlcard
2563 logical normalize,next
2565 double precision scal_bfac
2566 double precision xlink(4,0:4) /
2568 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2569 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2570 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2571 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2572 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2573 c print *, "WCHODZE"
2574 write (iout,*) "Calling read_dist_constr"
2575 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2577 restr_on_coord=.false.
2586 call card_concat(controlcard)
2587 next = index(controlcard,"NEXT").gt.0
2588 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2589 write (iout,*) "restr_type",restr_type
2590 call readi(controlcard,"NFRAG",nfrag_,0)
2591 call readi(controlcard,"NFRAG",nfrag_,0)
2592 call readi(controlcard,"NPAIR",npair_,0)
2593 call readi(controlcard,"NDIST",ndist_,0)
2594 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2595 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2596 if (restr_type.eq.10)
2597 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2598 if (restr_type.eq.12)
2599 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2600 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2601 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2602 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2603 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2604 normalize = index(controlcard,"NORMALIZE").gt.0
2605 write (iout,*) "WBOLTZD",wboltzd
2606 write (iout,*) "SCAL_PEAK",scal_peak
2607 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2608 write (iout,*) "IFRAG"
2610 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2612 write (iout,*) "IPAIR"
2614 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2616 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2618 & "Distance restraints as generated from reference structure"
2620 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2621 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2622 & ifrag_(2,i)=nstart_sup+nsup-1
2623 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2625 if (wfrag_(i).eq.0.0d0) cycle
2626 do j=ifrag_(1,i),ifrag_(2,i)-1
2627 do k=j+1,ifrag_(2,i)
2628 c write (iout,*) "j",j," k",k
2630 if (restr_type.eq.1) then
2636 forcon(nhpb)=wfrag_(i)
2637 else if (constr_dist.eq.2) then
2638 if (ddjk.le.dist_cut) then
2644 forcon(nhpb)=wfrag_(i)
2646 else if (restr_type.eq.3) then
2652 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2655 if (.not.out1file .or. me.eq.king)
2656 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2657 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2659 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2660 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2666 if (wpair_(i).eq.0.0d0) cycle
2674 do j=ifrag_(1,ii),ifrag_(2,ii)
2675 do k=ifrag_(1,jj),ifrag_(2,jj)
2677 if (restr_type.eq.1) then
2683 forcon(nhpb)=wpair_(i)
2684 else if (constr_dist.eq.2) then
2685 if (ddjk.le.dist_cut) then
2691 forcon(nhpb)=wpair_(i)
2693 else if (restr_type.eq.3) then
2699 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2702 if (.not.out1file .or. me.eq.king)
2703 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2704 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2706 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2707 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2714 write (iout,*) "Distance restraints as read from input"
2716 if (restr_type.eq.12) then
2717 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2718 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2719 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2720 & fordepth_peak(nhpb_peak+1),npeak
2721 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2722 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2723 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2724 c & fordepth_peak(nhpb_peak+1),npeak
2725 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2726 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2727 nhpb_peak=nhpb_peak+1
2728 irestr_type_peak(nhpb_peak)=12
2729 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2732 if (.not.out1file .or. me.eq.king)
2733 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2734 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2735 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2736 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2737 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2739 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2740 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2741 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2742 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2743 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2745 if (ibecarb_peak(nhpb_peak).eq.3) then
2746 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2747 else if (ibecarb_peak(nhpb_peak).eq.2) then
2748 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2749 else if (ibecarb_peak(nhpb_peak).eq.1) then
2750 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2751 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2753 else if (restr_type.eq.11) then
2754 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2755 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2756 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2757 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2759 irestr_type(nhpb)=11
2761 if (.not.out1file .or. me.eq.king)
2762 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2763 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2764 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2766 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2767 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2768 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2770 c if (ibecarb(nhpb).gt.0) then
2771 c ihpb(nhpb)=ihpb(nhpb)+nres
2772 c jhpb(nhpb)=jhpb(nhpb)+nres
2774 if (ibecarb(nhpb).eq.3) then
2775 ihpb(nhpb)=ihpb(nhpb)+nres
2776 else if (ibecarb(nhpb).eq.2) then
2777 ihpb(nhpb)=ihpb(nhpb)+nres
2778 else if (ibecarb(nhpb).eq.1) then
2779 ihpb(nhpb)=ihpb(nhpb)+nres
2780 jhpb(nhpb)=jhpb(nhpb)+nres
2782 else if (restr_type.eq.10) then
2783 c Cross-lonk Markov-like potential
2784 call card_concat(controlcard)
2785 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2786 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2788 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2789 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2790 if (index(controlcard,"ZL").gt.0) then
2792 else if (index(controlcard,"ADH").gt.0) then
2794 else if (index(controlcard,"PDH").gt.0) then
2796 else if (index(controlcard,"DSS").gt.0) then
2801 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2802 & xlink(1,link_type))
2803 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2804 & xlink(2,link_type))
2805 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2806 & xlink(3,link_type))
2807 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2808 & xlink(4,link_type))
2809 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2810 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2811 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2812 if (forcon(nhpb+1).le.0.0d0 .or.
2813 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2815 irestr_type(nhpb)=10
2816 if (ibecarb(nhpb).eq.3) then
2817 jhpb(nhpb)=jhpb(nhpb)+nres
2818 else if (ibecarb(nhpb).eq.2) then
2819 ihpb(nhpb)=ihpb(nhpb)+nres
2820 else if (ibecarb(nhpb).eq.1) then
2821 ihpb(nhpb)=ihpb(nhpb)+nres
2822 jhpb(nhpb)=jhpb(nhpb)+nres
2825 if (.not.out1file .or. me.eq.king)
2826 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2827 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2828 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2831 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2832 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2833 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2838 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2839 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2840 if (forcon(nhpb+1).gt.0.0d0) then
2842 if (dhpb1(nhpb).eq.0.0d0) then
2847 if (ibecarb(nhpb).eq.3) then
2848 jhpb(nhpb)=jhpb(nhpb)+nres
2849 else if (ibecarb(nhpb).eq.2) then
2850 ihpb(nhpb)=ihpb(nhpb)+nres
2851 else if (ibecarb(nhpb).eq.1) then
2852 ihpb(nhpb)=ihpb(nhpb)+nres
2853 jhpb(nhpb)=jhpb(nhpb)+nres
2855 if (dhpb(nhpb).eq.0.0d0)
2856 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2859 if (.not.out1file .or. me.eq.king)
2860 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2861 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2863 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2864 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2867 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2868 C if (forcon(nhpb+1).gt.0.0d0) then
2870 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2873 if (restr_type.eq.4) then
2874 write (iout,*) "The BFAC array"
2876 write (iout,'(i5,f10.5)') i,bfac(i)
2879 if (itype(i).eq.ntyp1) cycle
2881 if (itype(j).eq.ntyp1) cycle
2882 if (itype(i).eq.10) then
2887 if (itype(j).eq.10) then
2897 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2900 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2901 ihpb(nhpb)=i+nres*ii
2902 jhpb(nhpb)=j+nres*jj
2903 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2905 if (.not.out1file .or. me.eq.king) then
2906 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2907 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2908 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2912 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2913 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2914 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2924 if (restr_type.eq.5) then
2925 restr_on_coord=.true.
2927 if (itype(i).eq.ntyp1) cycle
2928 bfac(i)=(scal_bfac/bfac(i))**2
2937 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2938 & fordepthmax=fordepth(i)
2941 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2946 c-------------------------------------------------------------------------------
2947 subroutine read_constr_homology
2949 include 'DIMENSIONS'
2953 include 'COMMON.SETUP'
2954 include 'COMMON.CONTROL'
2955 include 'COMMON.HOMOLOGY'
2956 include 'COMMON.CHAIN'
2957 include 'COMMON.IOUNITS'
2959 include 'COMMON.QRESTR'
2960 include 'COMMON.GEO'
2961 include 'COMMON.INTERACT'
2962 include 'COMMON.NAMES'
2964 c For new homol impl
2966 include 'COMMON.VAR'
2969 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2971 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2972 c & sigma_odl_temp(maxres,maxres,max_template)
2974 character*24 model_ki_dist, model_ki_angle
2975 character*500 controlcard
2976 integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
2977 & ik,iistart,iishift
2982 c FP - Nov. 2014 Temporary specifications for new vars
2984 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2986 double precision, dimension (max_template,maxres) :: rescore
2987 double precision, dimension (max_template,maxres) :: rescore2
2988 double precision, dimension (max_template,maxres) :: rescore3
2989 double precision distal
2990 character*24 pdbfile,tpl_k_rescore
2991 c -----------------------------------------------------------------
2992 c Reading multiple PDB ref structures and calculation of retraints
2993 c not using pre-computed ones stored in files model_ki_{dist,angle}
2995 c -----------------------------------------------------------------
2998 c Alternative: reading from input
2999 call card_concat(controlcard)
3000 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3001 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3002 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3003 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3004 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3005 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3006 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3007 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3008 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3009 if(.not.read2sigma.and.start_from_model) then
3010 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3011 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3012 start_from_model=.false.
3014 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3015 & write(iout,*) 'START_FROM_MODELS is ON'
3016 if(start_from_model .and. rest) then
3017 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3018 write(iout,*) 'START_FROM_MODELS is OFF'
3019 write(iout,*) 'remove restart keyword from input'
3022 if (homol_nset.gt.1)then
3023 call card_concat(controlcard)
3024 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3025 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3026 write(iout,*) "iset homology_weight "
3028 write(iout,*) i,waga_homology(i)
3031 iset=mod(kolor,homol_nset)+1
3034 waga_homology(1)=1.0
3037 cd write (iout,*) "nnt",nnt," nct",nct
3044 c write(iout,*) 'nnt=',nnt,'nct=',nct
3047 do k=1,constr_homology
3060 if (read_homol_frag) then
3061 call read_klapaucjusz
3064 do k=1,constr_homology
3066 read(inp,'(a)') pdbfile
3067 if(me.eq.king .or. .not. out1file)
3068 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3069 & pdbfile(:ilen(pdbfile))
3070 open(ipdbin,file=pdbfile,status='old',err=33)
3072 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3073 & pdbfile(:ilen(pdbfile))
3076 c print *,'Begin reading pdb data'
3078 c Files containing res sim or local scores (former containing sigmas)
3081 write(kic2,'(bz,i2.2)') k
3083 tpl_k_rescore="template"//kic2//".sco"
3086 if (read2sigma) then
3087 call readpdb_template(k)
3092 c Distance restraints
3095 C Copy the coordinates from reference coordinates (?)
3099 c write (iout,*) "c(",j,i,") =",c(j,i)
3103 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3105 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3106 open (ientin,file=tpl_k_rescore,status='old')
3107 if (nnt.gt.1) rescore(k,1)=0.0d0
3108 do irec=nnt,nct ! loop for reading res sim
3109 if (read2sigma) then
3110 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3111 & rescore3_tmp,idomain_tmp
3113 idomain(k,i_tmp)=idomain_tmp
3114 rescore(k,i_tmp)=rescore_tmp
3115 rescore2(k,i_tmp)=rescore2_tmp
3116 rescore3(k,i_tmp)=rescore3_tmp
3117 if (.not. out1file .or. me.eq.king)
3118 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3119 & i_tmp,rescore2_tmp,rescore_tmp,
3120 & rescore3_tmp,idomain_tmp
3123 read (ientin,*,end=1401) rescore_tmp
3125 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3126 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3127 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3132 if (waga_dist.ne.0.0d0) then
3140 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3141 c write (iout,*) k,i,j,distal,dist2_cut
3143 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3144 & .and. distal.le.dist2_cut ) then
3150 c write (iout,*) "k",k
3151 c write (iout,*) "i",i," j",j," constr_homology",
3156 if (read2sigma) then
3159 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3161 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3162 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3163 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3165 if (odl(k,ii).le.dist_cut) then
3166 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3169 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3170 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3172 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3173 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3177 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3180 l_homo(k,ii)=.false.
3187 c Theta, dihedral and SC retraints
3189 if (waga_angle.gt.0.0d0) then
3190 c open (ientin,file=tpl_k_sigma_dih,status='old')
3191 c do irec=1,maxres-3 ! loop for reading sigma_dih
3192 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3193 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3194 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3195 c & sigma_dih(k,i+nnt-1)
3200 if (idomain(k,i).eq.0) then
3204 dih(k,i)=phiref(i) ! right?
3205 c read (ientin,*) sigma_dih(k,i) ! original variant
3206 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3207 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3208 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3209 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3210 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3212 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3213 & rescore(k,i-2)+rescore(k,i-3))/4.0
3214 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3215 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3216 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3217 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3218 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3219 c Instead of res sim other local measure of b/b str reliability possible
3220 if (sigma_dih(k,i).ne.0)
3221 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3222 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3227 if (waga_theta.gt.0.0d0) then
3228 c open (ientin,file=tpl_k_sigma_theta,status='old')
3229 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3230 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3231 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3232 c & sigma_theta(k,i+nnt-1)
3237 do i = nnt+2,nct ! right? without parallel.
3238 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3239 c do i=ithet_start,ithet_end ! with FG parallel.
3240 if (idomain(k,i).eq.0) then
3241 sigma_theta(k,i)=0.0
3244 thetatpl(k,i)=thetaref(i)
3245 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3246 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3247 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3248 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3249 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3250 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3251 & rescore(k,i-2))/3.0
3252 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3253 if (sigma_theta(k,i).ne.0)
3254 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3256 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3257 c rescore(k,i-2) ! right expression ?
3258 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3262 if (waga_d.gt.0.0d0) then
3263 c open (ientin,file=tpl_k_sigma_d,status='old')
3264 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3265 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3266 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3267 c & sigma_d(k,i+nnt-1)
3271 do i = nnt,nct ! right? without parallel.
3272 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3273 c do i=loc_start,loc_end ! with FG parallel.
3274 if (itype(i).eq.10) cycle
3275 if (idomain(k,i).eq.0 ) then
3282 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3283 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3284 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3285 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3286 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3287 if (sigma_d(k,i).ne.0)
3288 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3290 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3291 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3292 c read (ientin,*) sigma_d(k,i) ! 1st variant
3297 c remove distance restraints not used in any model from the list
3298 c shift data in all arrays
3300 if (waga_dist.ne.0.0d0) then
3306 if (ii_in_use(ii).eq.0.and.liiflag) then
3310 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3311 & .not.liiflag.and.ii.eq.lim_odl) then
3312 if (ii.eq.lim_odl) then
3313 iishift=ii-iistart+1
3318 do ki=iistart,lim_odl-iishift
3319 ires_homo(ki)=ires_homo(ki+iishift)
3320 jres_homo(ki)=jres_homo(ki+iishift)
3321 ii_in_use(ki)=ii_in_use(ki+iishift)
3322 do k=1,constr_homology
3323 odl(k,ki)=odl(k,ki+iishift)
3324 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3325 l_homo(k,ki)=l_homo(k,ki+iishift)
3329 lim_odl=lim_odl-iishift
3335 endif ! .not. klapaucjusz
3337 if (constr_homology.gt.0) call homology_partition
3338 if (constr_homology.gt.0) call init_int_table
3339 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3340 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3344 if (.not.out_template_restr) return
3345 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3346 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3347 write (iout,*) "Distance restraints from templates"
3349 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3350 & ii,ires_homo(ii),jres_homo(ii),
3351 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3352 & ki=1,constr_homology)
3354 write (iout,*) "Dihedral angle restraints from templates"
3356 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3357 & (rad2deg*dih(ki,i),
3358 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3360 write (iout,*) "Virtual-bond angle restraints from templates"
3362 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3363 & (rad2deg*thetatpl(ki,i),
3364 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3366 write (iout,*) "SC restraints from templates"
3368 write(iout,'(i5,100(4f8.2,4x))') i,
3369 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3370 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3373 c -----------------------------------------------------------------
3376 c----------------------------------------------------------------------
3378 subroutine flush(iu)
3383 subroutine flush(iu)
3388 c------------------------------------------------------------------------------
3389 subroutine copy_to_tmp(source)
3391 include "DIMENSIONS"
3392 include "COMMON.IOUNITS"
3393 character*(*) source
3394 character* 256 tmpfile
3398 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3399 inquire(file=tmpfile,exist=ex)
3401 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3402 & " to temporary directory..."
3403 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3404 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3408 c------------------------------------------------------------------------------
3409 subroutine move_from_tmp(source)
3411 include "DIMENSIONS"
3412 include "COMMON.IOUNITS"
3413 character*(*) source
3416 write (*,*) "Moving ",source(:ilen(source)),
3417 & " from temporary directory to working directory"
3418 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3419 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3422 c------------------------------------------------------------------------------
3423 subroutine random_init(seed)
3425 C Initialize random number generator
3428 include 'DIMENSIONS'
3431 logical OKRandom, prng_restart
3433 integer iseed_array(4)
3434 integer error_msg,ierr
3436 include 'COMMON.IOUNITS'
3437 include 'COMMON.TIME1'
3438 include 'COMMON.THREAD'
3439 include 'COMMON.SBRIDGE'
3440 include 'COMMON.CONTROL'
3441 include 'COMMON.MCM'
3442 include 'COMMON.MAP'
3443 include 'COMMON.HEADER'
3444 include 'COMMON.CSA'
3445 include 'COMMON.CHAIN'
3446 include 'COMMON.MUCA'
3448 include 'COMMON.FFIELD'
3449 include 'COMMON.SETUP'
3451 double precision seed,ran_number
3452 iseed=-dint(dabs(seed))
3453 if (iseed.eq.0) then
3454 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3455 & 'Random seed undefined. The program will stop.'
3456 write (*,'(/80(1h*)/20x,a/80(1h*))')
3457 & 'Random seed undefined. The program will stop.'
3459 call mpi_finalize(mpi_comm_world,ierr)
3461 stop 'Bad random seed.'
3464 if (fg_rank.eq.0) then
3468 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3469 OKRandom = prng_restart(me,iseed)
3472 tmp=65536.0d0**(4-i)
3473 iseed_array(i) = dint(seed/tmp)
3474 seed=seed-iseed_array(i)*tmp
3477 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3478 & (iseed_array(i),i=1,4)
3479 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3480 & (iseed_array(i),i=1,4)
3481 OKRandom = prng_restart(me,iseed_array)
3484 c r1 = prng_next(me)
3485 r1=ran_number(0.0D0,1.0D0)
3487 & write (iout,*) 'ran_num',r1
3488 if (r1.lt.0.0d0) OKRandom=.false.
3490 if (.not.OKRandom) then
3491 write (iout,*) 'PRNG IS NOT WORKING!!!'
3492 print *,'PRNG IS NOT WORKING!!!'
3495 call mpi_abort(mpi_comm_world,error_msg,ierr)
3498 write (iout,*) 'too many processors for parallel prng'
3499 write (*,*) 'too many processors for parallel prng'
3507 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3511 c----------------------------------------------------------------------
3512 subroutine read_klapaucjusz
3514 include 'DIMENSIONS'
3518 include 'COMMON.SETUP'
3519 include 'COMMON.CONTROL'
3520 include 'COMMON.HOMOLOGY'
3521 include 'COMMON.CHAIN'
3522 include 'COMMON.IOUNITS'
3524 include 'COMMON.GEO'
3525 include 'COMMON.INTERACT'
3526 include 'COMMON.NAMES'
3527 character*256 fragfile
3528 integer ninclust(maxclust),inclust(max_template,maxclust),
3529 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
3532 character*24 model_ki_dist, model_ki_angle
3533 character*500 controlcard
3534 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3535 & ik,ll,ii,kk,iistart,iishift,lim_xx
3536 double precision distal
3537 logical lprn /.true./
3543 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3544 double precision, dimension (max_template,maxres) :: rescore
3545 double precision, dimension (max_template,maxres) :: rescore2
3546 character*24 pdbfile,tpl_k_rescore
3549 c For new homol impl
3551 include 'COMMON.VAR'
3553 call getenv("FRAGFILE",fragfile)
3554 open(ientin,file=fragfile,status="old",err=10)
3555 read(ientin,*) constr_homology,nclust
3561 do k=1,constr_homology
3562 read(ientin,'(a)') pdbfile
3563 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3564 & pdbfile(:ilen(pdbfile))
3565 open(ipdbin,file=pdbfile,status='old',err=33)
3567 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3568 & pdbfile(:ilen(pdbfile))
3572 call readpdb_template(k)
3580 read(ientin,*) ninclust(i),nresclust(i)
3581 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3582 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3585 c Loop over clusters
3588 do ll = 1,ninclust(l)
3596 idomain(k,iresclust(i,l)+1) = 1
3598 idomain(k,iresclust(i,l)) = 1
3602 c Distance restraints
3605 C Copy the coordinates from reference coordinates (?)
3609 c write (iout,*) "c(",j,i,") =",c(j,i)
3612 call int_from_cart(.true.,.false.)
3613 call sc_loc_geom(.false.)
3615 thetaref(i)=theta(i)
3618 if (waga_dist.ne.0.0d0) then
3626 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3627 c write (iout,*) k,i,j,distal,dist2_cut
3629 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3630 & .and. distal.le.dist2_cut ) then
3636 c write (iout,*) "k",k
3637 c write (iout,*) "i",i," j",j," constr_homology",
3642 if (read2sigma) then
3645 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3647 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3648 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3649 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3651 if (odl(k,ii).le.dist_cut) then
3652 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3655 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3656 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3658 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3659 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3663 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3666 c l_homo(k,ii)=.false.
3673 c Theta, dihedral and SC retraints
3675 if (waga_angle.gt.0.0d0) then
3677 if (idomain(k,i).eq.0) then
3678 c sigma_dih(k,i)=0.0
3682 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3683 & rescore(k,i-2)+rescore(k,i-3))/4.0
3684 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3685 c & " sigma_dihed",sigma_dih(k,i)
3686 if (sigma_dih(k,i).ne.0)
3687 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3692 if (waga_theta.gt.0.0d0) then
3694 if (idomain(k,i).eq.0) then
3695 c sigma_theta(k,i)=0.0
3698 thetatpl(k,i)=thetaref(i)
3699 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3700 & rescore(k,i-2))/3.0
3701 if (sigma_theta(k,i).ne.0)
3702 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3706 if (waga_d.gt.0.0d0) then
3708 if (itype(i).eq.10) cycle
3709 if (idomain(k,i).eq.0 ) then
3716 sigma_d(k,i)=rescore(k,i)
3717 if (sigma_d(k,i).ne.0)
3718 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3719 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3725 c remove distance restraints not used in any model from the list
3726 c shift data in all arrays
3728 if (waga_dist.ne.0.0d0) then
3734 if (ii_in_use(ii).eq.0.and.liiflag) then
3738 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3739 & .not.liiflag.and.ii.eq.lim_odl) then
3740 if (ii.eq.lim_odl) then
3741 iishift=ii-iistart+1
3746 do ki=iistart,lim_odl-iishift
3747 ires_homo(ki)=ires_homo(ki+iishift)
3748 jres_homo(ki)=jres_homo(ki+iishift)
3749 ii_in_use(ki)=ii_in_use(ki+iishift)
3750 do k=1,constr_homology
3751 odl(k,ki)=odl(k,ki+iishift)
3752 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3753 l_homo(k,ki)=l_homo(k,ki+iishift)
3757 lim_odl=lim_odl-iishift
3764 10 stop "Error in fragment file"