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 readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
175 call reada(controlcard,"T_BATH",t_bath,300.0d0)
176 minim=(index(controlcard,'MINIMIZE').gt.0)
177 dccart=(index(controlcard,'CART').gt.0)
178 overlapsc=(index(controlcard,'OVERLAP').gt.0)
179 overlapsc=.not.overlapsc
180 searchsc=(index(controlcard,'SEARCHSC').gt.0 .and.
181 & index(controlcard,'NOSEARCHSC').eq.0)
182 sideadd=(index(controlcard,'SIDEADD').gt.0)
183 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
184 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
185 outpdb=(index(controlcard,'PDBOUT').gt.0)
186 outmol2=(index(controlcard,'MOL2OUT').gt.0)
187 pdbref=(index(controlcard,'PDBREF').gt.0)
188 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
189 indpdb=index(controlcard,'PDBSTART')
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 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
299 extconf=(index(controlcard,'EXTCONF').gt.0)
300 if (start_from_model) then
304 i2ndstr=index(controlcard,'USE_SEC_PRED')
305 gradout=index(controlcard,'GRADOUT').gt.0
306 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
307 C DISTCHAINMAX become obsolete for periodic boundry condition
308 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
309 C Reading the dimensions of box in x,y,z coordinates
310 call reada(controlcard,'BOXX',boxxsize,100.0d0)
311 call reada(controlcard,'BOXY',boxysize,100.0d0)
312 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
313 write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
314 c Cutoff range for interactions
315 call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
316 call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
317 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
318 write (iout,*) "Cutoff on interactions",r_cut_int
320 & "Cutoff in switching short and long range interactions in RESPA",
322 write (iout,*) "lambda in switch function",rlamb
323 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
324 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
325 if (lipthick.gt.0.0d0) then
326 bordliptop=(boxzsize+lipthick)/2.0
327 bordlipbot=bordliptop-lipthick
329 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
330 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
331 buflipbot=bordlipbot+lipbufthick
332 bufliptop=bordliptop-lipbufthick
333 if ((lipbufthick*2.0d0).gt.lipthick)
334 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
336 write(iout,*) "bordliptop=",bordliptop
337 write(iout,*) "bordlipbot=",bordlipbot
338 write(iout,*) "bufliptop=",bufliptop
339 write(iout,*) "buflipbot=",buflipbot
340 write (iout,*) "SHIELD MODE",shield_mode
341 if (TUBElog.gt.0) then
342 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
343 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
344 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
345 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
346 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
347 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
348 buftubebot=bordtubebot+tubebufthick
349 buftubetop=bordtubetop-tubebufthick
351 if (me.eq.king .or. .not.out1file )
352 & write (iout,*) "DISTCHAINMAX",distchainmax
354 if(me.eq.king.or..not.out1file)
355 & write (iout,'(2a)') diagmeth(kdiag),
356 & ' routine used to diagonalize matrices.'
359 c--------------------------------------------------------------------------
360 subroutine read_REMDpar
366 include 'COMMON.IOUNITS'
367 include 'COMMON.TIME1'
370 include 'COMMON.LANGEVIN'
373 include 'COMMON.LANGEVIN.lang0.5diag'
375 include 'COMMON.LANGEVIN.lang0'
378 include 'COMMON.INTERACT'
379 include 'COMMON.NAMES'
381 include 'COMMON.REMD'
382 include 'COMMON.CONTROL'
383 include 'COMMON.SETUP'
385 character*320 controlcard
386 character*3200 controlcard1
387 integer iremd_m_total,i
389 if(me.eq.king.or..not.out1file)
390 & write (iout,*) "REMD setup"
392 call card_concat(controlcard)
393 call readi(controlcard,"NREP",nrep,3)
394 call readi(controlcard,"NSTEX",nstex,1000)
395 call reada(controlcard,"RETMIN",retmin,10.0d0)
396 call reada(controlcard,"RETMAX",retmax,1000.0d0)
397 mremdsync=(index(controlcard,'SYNC').gt.0)
398 call readi(controlcard,"NSYN",i_sync_step,100)
399 restart1file=(index(controlcard,'REST1FILE').gt.0)
400 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
401 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
402 if(max_cache_traj_use.gt.max_cache_traj)
403 & max_cache_traj_use=max_cache_traj
404 if(me.eq.king.or..not.out1file) then
405 cd if (traj1file) then
406 crc caching is in testing - NTWX is not ignored
407 cd write (iout,*) "NTWX value is ignored"
408 cd write (iout,*) " trajectory is stored to one file by master"
409 cd write (iout,*) " before exchange at NSTEX intervals"
411 write (iout,*) "NREP= ",nrep
412 write (iout,*) "NSTEX= ",nstex
413 write (iout,*) "SYNC= ",mremdsync
414 write (iout,*) "NSYN= ",i_sync_step
415 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
418 if (index(controlcard,'TLIST').gt.0) then
420 call card_concat(controlcard1)
421 read(controlcard1,*) (remd_t(i),i=1,nrep)
422 if(me.eq.king.or..not.out1file)
423 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
426 if (index(controlcard,'MLIST').gt.0) then
428 call card_concat(controlcard1)
429 read(controlcard1,*) (remd_m(i),i=1,nrep)
430 if(me.eq.king.or..not.out1file) then
431 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
434 iremd_m_total=iremd_m_total+remd_m(i)
436 write (iout,*) 'Total number of replicas ',iremd_m_total
441 & "Adaptive (PMF-biased) umbrella sampling will be run"
444 if(me.eq.king.or..not.out1file)
445 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
448 c--------------------------------------------------------------------------
449 subroutine read_MDpar
455 include 'COMMON.IOUNITS'
456 include 'COMMON.TIME1'
458 include 'COMMON.QRESTR'
460 include 'COMMON.LANGEVIN'
463 include 'COMMON.LANGEVIN.lang0.5diag'
465 include 'COMMON.LANGEVIN.lang0'
468 include 'COMMON.INTERACT'
469 include 'COMMON.NAMES'
471 include 'COMMON.SETUP'
472 include 'COMMON.CONTROL'
473 include 'COMMON.SPLITELE'
474 include 'COMMON.FFIELD'
476 character*320 controlcard
480 call card_concat(controlcard)
481 call readi(controlcard,"NSTEP",n_timestep,1000000)
482 call readi(controlcard,"NTWE",ntwe,100)
483 call readi(controlcard,"NTWX",ntwx,1000)
484 call reada(controlcard,"DT",d_time,1.0d-1)
485 call reada(controlcard,"DVMAX",dvmax,2.0d1)
486 call reada(controlcard,"DAMAX",damax,1.0d1)
487 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
488 call readi(controlcard,"LANG",lang,0)
489 RESPA = index(controlcard,"RESPA") .gt. 0
490 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
491 ntime_split0=ntime_split
492 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
493 ntime_split0=ntime_split
494 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
495 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
496 rest = index(controlcard,"REST").gt.0
497 tbf = index(controlcard,"TBF").gt.0
498 usampl = index(controlcard,"USAMPL").gt.0
499 scale_umb = index(controlcard,"SCALE_UMB").gt.0
500 adaptive = index(controlcard,"ADAPTIVE").gt.0
501 mdpdb = index(controlcard,"MDPDB").gt.0
502 call reada(controlcard,"T_BATH",t_bath,300.0d0)
503 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
504 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
505 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
506 if (count_reset_moment.eq.0) count_reset_moment=1000000000
507 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
508 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
509 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
510 if (count_reset_vel.eq.0) count_reset_vel=1000000000
511 large = index(controlcard,"LARGE").gt.0
512 print_compon = index(controlcard,"PRINT_COMPON").gt.0
513 rattle = index(controlcard,"RATTLE").gt.0
514 preminim = index(controlcard,"PREMINIM").gt.0
515 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
519 write (iout,*) "PREMINIM ",preminim
521 if (index(controlcard,'CART').gt.0) dccart=.true.
522 write (iout,*) "dccart ",dccart
523 write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
524 if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
526 c if performing umbrella sampling, fragments constrained are read from the fragment file
529 write (iout,*) "Umbrella sampling will be run"
530 if (scale_umb.and.adaptive) then
531 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
532 write (iout,*) "Select one of those and re-run the job."
535 if (scale_umb) write (iout,*)
536 &"Umbrella-restraint force constants will be scaled by temperature"
540 if(me.eq.king.or..not.out1file) then
542 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
544 write (iout,'(a)') "The units are:"
545 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
546 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
547 & " acceleration: angstrom/(48.9 fs)**2"
548 write (iout,'(a)') "energy: kcal/mol, temperature: K"
550 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
551 write (iout,'(a60,f10.5,a)')
552 & "Initial time step of numerical integration:",d_time,
554 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
555 write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
557 write(iout,'(a60,i5)')"Frequency of updating interaction list",
560 write (iout,'(2a,i4,a)')
561 & "A-MTS algorithm used; initial time step for fast-varying",
562 & " short-range forces split into",ntime_split," steps."
563 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
564 & r_cut_respa," lambda",rlamb
566 write (iout,'(2a,f10.5)')
567 & "Maximum acceleration threshold to reduce the time step",
568 & "/increase split number:",damax
569 write (iout,'(2a,f10.5)')
570 & "Maximum predicted energy drift to reduce the timestep",
571 & "/increase split number:",edriftmax
572 write (iout,'(a60,f10.5)')
573 & "Maximum velocity threshold to reduce velocities:",dvmax
574 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
575 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
576 if (rattle) write (iout,'(a60)')
577 & "Rattle algorithm used to constrain the virtual bonds"
581 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
582 call reada(controlcard,"RWAT",rwat,1.4d0)
583 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
584 surfarea=index(controlcard,"SURFAREA").gt.0
585 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
586 if(me.eq.king.or..not.out1file)then
587 write (iout,'(/a,$)') "Langevin dynamics calculation"
590 & " with direct integration of Langevin equations"
591 else if (lang.eq.2) then
592 write (iout,'(a/)') " with TINKER stochasic MD integrator"
593 else if (lang.eq.3) then
594 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
595 else if (lang.eq.4) then
596 write (iout,'(a/)') " in overdamped mode"
598 write (iout,'(//a,i5)')
599 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
602 write (iout,'(a60,f10.5)') "Temperature:",t_bath
603 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
604 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
605 write (iout,'(a60,f10.5)')
606 & "Scaling factor of the friction forces:",scal_fric
607 if (surfarea) write (iout,'(2a,i10,a)')
608 & "Friction coefficients will be scaled by solvent-accessible",
609 & " surface area every",reset_fricmat," steps."
611 c Calculate friction coefficients and bounds of stochastic forces
612 eta=6*pi*cPoise*etawat
613 if(me.eq.king.or..not.out1file)
614 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
616 gamp=scal_fric*(pstok+rwat)*eta
617 stdfp=dsqrt(2*Rb*t_bath/d_time)
619 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
620 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
622 if(me.eq.king.or..not.out1file)then
623 write (iout,'(/2a/)')
624 & "Radii of site types and friction coefficients and std's of",
625 & " stochastic forces of fully exposed sites"
626 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
628 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
629 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
633 if(me.eq.king.or..not.out1file)then
634 write (iout,'(a)') "Berendsen bath calculation"
635 write (iout,'(a60,f10.5)') "Temperature:",t_bath
636 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
638 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
639 & count_reset_moment," steps"
641 & write (iout,'(a,i10,a)')
642 & "Velocities will be reset at random every",count_reset_vel,
646 if(me.eq.king.or..not.out1file)
647 & write (iout,'(a31)') "Microcanonical mode calculation"
649 if(me.eq.king.or..not.out1file)then
650 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
652 write(iout,*) "MD running with constraints."
653 write(iout,*) "Equilibration time ", eq_time, " mtus."
654 write(iout,*) "Constraining ", nfrag," fragments."
655 write(iout,*) "Length of each fragment, weight and q0:"
657 write (iout,*) "Set of restraints #",iset
659 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
660 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
662 write(iout,*) "constraints between ", npair, "fragments."
663 write(iout,*) "constraint pairs, weights and q0:"
665 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
666 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
668 write(iout,*) "angle constraints within ", nfrag_back,
669 & "backbone fragments."
671 write(iout,*) "fragment, weights, q0:"
673 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
674 & ifrag_back(2,i,iset),
675 & wfrag_back(1,i,iset),qin_back(1,i,iset),
676 & wfrag_back(2,i,iset),qin_back(2,i,iset),
677 & wfrag_back(3,i,iset),qin_back(3,i,iset)
680 write(iout,*) "fragment, weights:"
682 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
683 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
684 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
688 iset=mod(kolor,nset)+1
691 if(me.eq.king.or..not.out1file)
692 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
695 c------------------------------------------------------------------------------
698 C Read molecular data.
704 integer error_msg,ierror,ierr,ierrcode
706 include 'COMMON.IOUNITS'
709 include 'COMMON.INTERACT'
710 include 'COMMON.LOCAL'
711 include 'COMMON.NAMES'
712 include 'COMMON.CHAIN'
713 include 'COMMON.FFIELD'
714 include 'COMMON.SBRIDGE'
715 include 'COMMON.HEADER'
716 include 'COMMON.CONTROL'
717 include 'COMMON.SAXS'
718 include 'COMMON.DBASE'
719 include 'COMMON.THREAD'
720 include 'COMMON.CONTACTS'
721 include 'COMMON.TORCNSTR'
722 include 'COMMON.TIME1'
723 include 'COMMON.BOUNDS'
725 include 'COMMON.SETUP'
726 include 'COMMON.SHIELD'
727 character*4 sequence(maxres)
729 double precision x(maxvar)
730 character*256 pdbfile
731 character*400 weightcard
732 character*80 weightcard_t,ucase
733 integer itype_pdb(maxres)
734 common /pizda/ itype_pdb
735 logical seq_comp,fail
736 double precision energia(0:n_ene)
737 double precision secprob(3,maxdih_constr)
739 double precision phihel,phibet,sigmahel,sigmabet
740 integer iti,nsi,maxsi
744 integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
745 double precision sumv
747 C Read PDB structure if applicable
749 if (indpdb.gt.0 .or. pdbref) then
750 read(inp,'(a)') pdbfile
751 if(me.eq.king.or..not.out1file)
752 & write (iout,'(2a)') 'PDB data will be read from file ',
753 & pdbfile(:ilen(pdbfile))
754 open(ipdbin,file=pdbfile,status='old',err=33)
756 33 write (iout,'(a)') 'Error opening PDB file.'
767 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
768 & (crefjlee(j,i+nres),j=1,3)
771 if(me.eq.king.or..not.out1file)
772 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
773 & ' nstart_sup=',nstart_sup
775 itype_pdb(i)=itype(i)
779 nct=nstart_sup+nsup-1
780 call contact(.false.,ncont_ref,icont_ref,co)
783 if(me.eq.king.or..not.out1file)
784 & write(iout,*)'Adding sidechains'
788 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
791 do while (fail.and.nsi.le.maxsi)
792 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
795 if(fail) write(iout,*)'Adding sidechain failed for res ',
796 & i,' after ',nsi,' trials'
801 if (indpdb.eq.0) then
802 C Read sequence if not taken from the pdb file.
804 c print *,'nres=',nres
805 if (iscode.gt.0) then
806 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
808 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
810 C Convert sequence to numeric code
812 itype(i)=rescode(i,sequence(i),iscode)
816 c print '(20i4)',(itype(i),i=1,nres)
819 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
821 if (itype(i).eq.ntyp1) then
825 else if (iabs(itype(i+1)).ne.20) then
827 else if (iabs(itype(i)).ne.20) then
834 if(me.eq.king.or..not.out1file)then
835 write (iout,*) "ITEL"
837 write (iout,*) i,itype(i),itel(i)
839 c print *,'Call Read_Bridge.'
843 cd print *,'NNT=',NNT,' NCT=',NCT
844 call seq2chains(nres,itype,nchain,chain_length,chain_border,
847 chain_border1(2,1)=chain_border(2,1)+1
849 chain_border1(1,i)=chain_border(1,i)-1
850 chain_border1(2,i)=chain_border(2,i)+1
852 if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
853 chain_border1(2,nchain)=nres
854 write(iout,*) "nres",nres," nchain",nchain
856 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
857 & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
859 call chain_symmetry(nchain,nres,itype,chain_border,
860 & chain_length,npermchain,tabpermchain)
862 c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
865 write(iout,*) "residue permutations"
867 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
870 if (itype(1).eq.ntyp1) nnt=2
871 if (itype(nres).eq.ntyp1) nct=nct-1
872 write (iout,*) "nnt",nnt," nct",nct
875 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
876 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
878 c print*, 'init_dfa_vars finished!'
880 c print*, 'read_dfa_info finished!'
884 if(me.eq.king.or..not.out1file)
885 & write (iout,'(a,i3)') 'nsup=',nsup
887 if (nsup.le.(nct-nnt+1)) then
888 do i=0,nct-nnt+1-nsup
889 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
895 & 'Error - sequences to be superposed do not match.'
898 do i=0,nsup-(nct-nnt+1)
899 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
901 nstart_sup=nstart_sup+i
907 & 'Error - sequences to be superposed do not match.'
910 if (nsup.eq.0) nsup=nct-nnt
911 if (nstart_sup.eq.0) nstart_sup=nnt
912 if (nstart_seq.eq.0) nstart_seq=nnt
913 if(me.eq.king.or..not.out1file)
914 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
915 & ' nstart_seq=',nstart_seq
918 C 8/13/98 Set limits to generating the dihedral angles
923 read (inp,*) ndih_constr
924 write (iout,*) "ndih_constr",ndih_constr
925 if (ndih_constr.gt.0) then
928 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
931 if(me.eq.king.or..not.out1file)then
934 & 'There are',ndih_constr,' restraints on gamma angles.'
936 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
943 phi0(i)=deg2rad*phi0(i)
944 drange(i)=deg2rad*drange(i)
946 C if(me.eq.king.or..not.out1file)
947 C & write (iout,*) 'FTORS',ftors
950 phibound(1,ii) = phi0(i)-drange(i)
951 phibound(2,ii) = phi0(i)+drange(i)
954 if (me.eq.king .or. .not.out1file) then
956 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
958 write (iout,'(a3,i5,2f10.1)')
959 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
964 else if (ndih_constr.lt.0) then
966 call card_concat(weightcard)
967 call reada(weightcard,"PHIHEL",phihel,50.0D0)
968 call reada(weightcard,"PHIBET",phibet,180.0D0)
969 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
970 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
971 call reada(weightcard,"WDIHC",wdihc,0.591D0)
972 write (iout,*) "Weight of dihedral angle restraints",wdihc
973 read(inp,'(9x,3f7.3)')
974 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
975 write (iout,*) "The secprob array"
977 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
981 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
982 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
983 ndih_constr=ndih_constr+1
984 idih_constr(ndih_constr)=i
987 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
988 sumv=sumv+vpsipred(j,ndih_constr)
991 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
993 phibound(1,ndih_constr)=phihel*deg2rad
994 phibound(2,ndih_constr)=phibet*deg2rad
995 sdihed(1,ndih_constr)=sigmahel*deg2rad
996 sdihed(2,ndih_constr)=sigmabet*deg2rad
1000 if(me.eq.king.or..not.out1file)then
1003 & 'There are',ndih_constr,
1004 & ' bimodal restraints on gamma angles.'
1006 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
1007 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
1008 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
1009 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1010 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1011 & (vpsipred(j,i),j=1,3)
1018 C first setting the theta boundaries to 0 to pi
1019 C this mean that there is no energy penalty for any angle occuring this can be applied
1020 C for generate random conformation but is not implemented in this way
1023 C thetabound(2,i)=pi
1025 C begin reading theta constrains this is quartic constrains allowing to
1026 C have smooth second derivative
1027 if (with_theta_constr) then
1028 C with_theta_constr is keyword allowing for occurance of theta constrains
1029 read (inp,*) ntheta_constr
1030 C ntheta_constr is the number of theta constrains
1031 if (ntheta_constr.gt.0) then
1032 C read (inp,*) ftors
1033 read (inp,*) (itheta_constr(i),theta_constr0(i),
1034 & theta_drange(i),for_thet_constr(i),
1035 & i=1,ntheta_constr)
1036 C the above code reads from 1 to ntheta_constr
1037 C itheta_constr(i) residue i for which is theta_constr
1038 C theta_constr0 the global minimum value
1039 C theta_drange is range for which there is no energy penalty
1040 C for_thet_constr is the force constant for quartic energy penalty
1042 if(me.eq.king.or..not.out1file)then
1044 & 'There are',ntheta_constr,' constraints on phi angles.'
1045 do i=1,ntheta_constr
1046 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1048 & for_thet_constr(i)
1051 do i=1,ntheta_constr
1052 theta_constr0(i)=deg2rad*theta_constr0(i)
1053 theta_drange(i)=deg2rad*theta_drange(i)
1055 C if(me.eq.king.or..not.out1file)
1056 C & write (iout,*) 'FTORS',ftors
1057 C do i=1,ntheta_constr
1058 C ii = itheta_constr(i)
1059 C thetabound(1,ii) = phi0(i)-drange(i)
1060 C thetabound(2,ii) = phi0(i)+drange(i)
1062 endif ! ntheta_constr.gt.0
1063 endif! with_theta_constr
1065 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1066 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1067 c--- Zscore rms -------
1068 if (nz_start.eq.0) nz_start=nnt
1069 if (nz_end.eq.0 .and. nsup.gt.0) then
1071 else if (nz_end.eq.0) then
1074 if(me.eq.king.or..not.out1file)then
1075 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1076 write (iout,*) 'IZ_SC=',iz_sc
1078 c----------------------
1081 if (.not.pdbref) then
1082 call read_angles(inp,*38)
1085 38 write (iout,'(a)') 'Error reading reference structure.'
1087 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1088 stop 'Error reading reference structure'
1090 39 call chainbuild_extconf
1092 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1101 call contact(.true.,ncont_ref,icont_ref,co)
1105 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1107 if (constr_dist.gt.0) call read_dist_constr
1108 c write (iout,*) "After read_dist_constr nhpb",nhpb
1109 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1111 call NMRpeak_partition
1112 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1114 if(me.eq.king.or..not.out1file)
1115 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1118 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1120 if(me.eq.king.or..not.out1file)
1121 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1122 & icont_ref(1,i),' ',
1123 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1126 write (iout,*) "calling read_saxs_consrtr",nsaxs
1127 if (nsaxs.gt.0) call read_saxs_constr
1129 if (constr_homology.gt.0) then
1130 call read_constr_homology
1131 if (indpdb.gt.0 .or. pdbref) then
1134 c(j,i)=crefjlee(j,i)
1135 cref(j,i)=crefjlee(j,i)
1140 write (iout,*) "sc_loc_geom: Array C"
1142 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1143 & (c(j,i+nres),j=1,3)
1145 write (iout,*) "Array Cref"
1147 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1148 & (cref(j,i+nres),j=1,3)
1151 call int_from_cart1(.false.)
1152 call sc_loc_geom(.false.)
1154 thetaref(i)=theta(i)
1159 dc(j,i)=c(j,i+1)-c(j,i)
1160 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1165 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1166 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1171 if (start_from_model) then
1174 read(inp,'(a)',end=332,err=332) pdbfile
1175 if (me.eq.king .or. .not. out1file)
1176 & write (iout,'(a,5x,a)') 'Opening PDB file',
1177 & pdbfile(:ilen(pdbfile))
1178 open(ipdbin,file=pdbfile,status='old',err=336)
1180 336 write (iout,'(a,5x,a)') 'Error opening PDB file',
1181 & pdbfile(:ilen(pdbfile))
1189 if (nres.ge.nres_temp) then
1190 nmodel_start=nmodel_start+1
1191 pdbfiles_chomo(nmodel_start)=pdbfile
1194 chomo(j,i,nmodel_start)=c(j,i)
1198 if (me.eq.king .or. .not. out1file)
1199 & write (iout,'(a,2i5,1x,a)')
1200 & "Different number of residues",nres_temp,nres,
1206 if (nmodel_start.eq.0) then
1207 if (me.eq.king .or. .not. out1file)
1208 & write (iout,'(a)')
1209 & "No valid starting model found START_FROM_MODELS is OFF"
1210 start_from_model=.false.
1212 write (iout,*) "nmodel_start",nmodel_start
1218 c write (iout,*) "iranconf",iranconf," extconf",extconf,
1219 c & " start_from_models",start_from_model
1220 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1221 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1222 & modecalc.ne.10) then
1223 C If input structure hasn't been supplied from the PDB file read or generate
1225 if (iranconf.eq.0 .and. .not. extconf .and. .not.
1226 & start_from_model) then
1227 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1228 & write (iout,'(a)') 'Initial geometry will be read in.'
1230 read(inp,'(8f10.5)',end=36,err=36)
1231 & ((c(l,k),l=1,3),k=1,nres),
1232 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1233 c write (iout,*) "Exit READ_CART"
1234 c write (iout,'(8f10.5)')
1235 c & ((c(l,k),l=1,3),k=1,nres),
1236 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1240 dc(j,i)=c(j,i+1)-c(j,i)
1241 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1245 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1247 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1248 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1253 c dc_norm(j,i+nres)=0.0d0
1257 call int_from_cart1(.true.)
1258 write (iout,*) "Finish INT_TO_CART"
1259 c write (iout,*) "DC"
1261 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1262 c & (dc(j,i+nres),j=1,3)
1268 call read_angles(inp,*36)
1270 call chainbuild_extconf
1273 36 write (iout,'(a)') 'Error reading angle file.'
1275 call mpi_finalize( MPI_COMM_WORLD,IERR )
1277 stop 'Error reading angle file.'
1279 else if (extconf) then
1280 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1281 & write (iout,'(a)') 'Extended chain initial geometry.'
1283 theta(i)=90d0*deg2rad
1286 phi(i)=180d0*deg2rad
1289 alph(i)=110d0*deg2rad
1292 omeg(i)=-120d0*deg2rad
1293 if (itype(i).le.0) omeg(i)=-omeg(i)
1296 call chainbuild_extconf
1298 if(me.eq.king.or..not.out1file)
1299 & write (iout,'(a)') 'Random-generated initial geometry.'
1302 if (me.eq.king .or. fg_rank.eq.0 .and. (
1303 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1307 call gen_rand_conf(itmp,*30)
1309 30 write (iout,*) 'Failed to generate random conformation',
1310 & ', itrial=',itrial
1311 write (*,*) 'Processor:',me,
1312 & ' Failed to generate random conformation',
1322 write (iout,'(a,i3,a)') 'Processor:',me,
1323 & ' error in generating random conformation.'
1324 write (*,'(a,i3,a)') 'Processor:',me,
1325 & ' error in generating random conformation.'
1328 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1333 & ' error in generating random conformation.'
1338 elseif (modecalc.eq.4) then
1339 read (inp,'(a)') intinname
1340 open (intin,file=intinname,status='old',err=333)
1341 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1342 & write (iout,'(a)') 'intinname',intinname
1343 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1345 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1347 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1349 stop 'Error opening angle file.'
1353 C Generate distance constraints, if the PDB structure is to be regularized.
1354 if (nthread.gt.0) then
1355 call read_threadbase
1358 if (me.eq.king .or. .not. out1file)
1360 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1361 write (iout,'(/a,i3,a)')
1362 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1363 write (iout,'(20i4)') (iss(i),i=1,ns)
1365 write(iout,*)"Running with dynamic disulfide-bond formation"
1367 write (iout,'(/a/)') 'Pre-formed links are:'
1373 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1374 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1380 if (ns.gt.0.and.dyn_ss) then
1384 forcon(i-nss)=forcon(i)
1391 dyn_ss_mask(iss(i))=.true.
1396 c write (iout,*) "DC"
1398 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1399 c & (dc(j,i+nres),j=1,3)
1401 if (i2ndstr.gt.0) call secstrp2dihc
1402 c call geom_to_var(nvar,x)
1403 c call etotal(energia(0))
1404 c call enerprint(energia(0))
1405 c call briefout(0,etot)
1407 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1408 cd write (iout,'(a)') 'Variable list:'
1409 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1411 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1412 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1413 & 'Processor',myrank,': end reading molecular data.'
1418 c--------------------------------------------------------------------------
1419 logical function seq_comp(itypea,itypeb,length)
1421 integer length,itypea(length),itypeb(length)
1424 if (itypea(i).ne.itypeb(i)) then
1432 c-----------------------------------------------------------------------------
1433 subroutine read_bridge
1434 C Read information about disulfide bridges.
1436 include 'DIMENSIONS'
1441 include 'COMMON.IOUNITS'
1442 include 'COMMON.GEO'
1443 include 'COMMON.VAR'
1444 include 'COMMON.INTERACT'
1445 include 'COMMON.LOCAL'
1446 include 'COMMON.NAMES'
1447 include 'COMMON.CHAIN'
1448 include 'COMMON.FFIELD'
1449 include 'COMMON.SBRIDGE'
1450 include 'COMMON.HEADER'
1451 include 'COMMON.CONTROL'
1452 include 'COMMON.DBASE'
1453 include 'COMMON.THREAD'
1454 include 'COMMON.TIME1'
1455 include 'COMMON.SETUP'
1457 C Read bridging residues.
1458 read (inp,*) ns,(iss(i),i=1,ns)
1459 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
1466 if(me.eq.king.or..not.out1file)
1467 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1468 C Check whether the specified bridging residues are cystines.
1470 if (itype(iss(i)).ne.1) then
1471 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1472 & 'Do you REALLY think that the residue ',
1473 & restyp(itype(iss(i))),i,
1474 & ' can form a disulfide bridge?!!!'
1475 write (*,'(2a,i3,a)')
1476 & 'Do you REALLY think that the residue ',
1477 & restyp(itype(iss(i))),i,
1478 & ' can form a disulfide bridge?!!!'
1480 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1485 C Read preformed bridges.
1487 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1489 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1492 C Check if the residues involved in bridges are in the specified list of
1493 C bridging residues.
1496 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1497 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1498 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1499 & ' contains residues present in other pairs.'
1500 write (*,'(a,i3,a)') 'Disulfide pair',i,
1501 & ' contains residues present in other pairs.'
1503 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1509 if (ihpb(i).eq.iss(j)) goto 10
1511 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1514 if (jhpb(i).eq.iss(j)) goto 20
1516 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1522 ihpb(i)=ihpb(i)+nres
1523 jhpb(i)=jhpb(i)+nres
1529 c----------------------------------------------------------------------------
1530 subroutine read_x(kanal,*)
1532 include 'DIMENSIONS'
1533 include 'COMMON.GEO'
1534 include 'COMMON.VAR'
1535 include 'COMMON.CHAIN'
1536 include 'COMMON.IOUNITS'
1537 include 'COMMON.CONTROL'
1538 include 'COMMON.LOCAL'
1539 include 'COMMON.INTERACT'
1540 integer i,j,k,l,kanal
1541 c Read coordinates from input
1543 read(kanal,'(8f10.5)',end=10,err=10)
1544 & ((c(l,k),l=1,3),k=1,nres),
1545 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1548 c(j,2*nres)=c(j,nres)
1550 call int_from_cart1(.false.)
1553 dc(j,i)=c(j,i+1)-c(j,i)
1554 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1558 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1560 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1561 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1569 c----------------------------------------------------------------------------
1570 subroutine read_threadbase
1572 include 'DIMENSIONS'
1573 include 'COMMON.IOUNITS'
1574 include 'COMMON.GEO'
1575 include 'COMMON.VAR'
1576 include 'COMMON.INTERACT'
1577 include 'COMMON.LOCAL'
1578 include 'COMMON.NAMES'
1579 include 'COMMON.CHAIN'
1580 include 'COMMON.FFIELD'
1581 include 'COMMON.SBRIDGE'
1582 include 'COMMON.HEADER'
1583 include 'COMMON.CONTROL'
1584 include 'COMMON.DBASE'
1585 include 'COMMON.THREAD'
1586 include 'COMMON.TIME1'
1588 double precision dist
1589 C Read pattern database for threading.
1590 read (icbase,*) nseq
1592 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1593 & nres_base(2,i),nres_base(3,i)
1594 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1596 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1597 c & nres_base(2,i),nres_base(3,i)
1598 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1602 if (weidis.eq.0.0D0) weidis=0.1D0
1611 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1612 write (iout,'(a,i5)') 'nexcl: ',nexcl
1613 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1616 c------------------------------------------------------------------------------
1617 subroutine setup_var
1619 include 'DIMENSIONS'
1620 include 'COMMON.IOUNITS'
1621 include 'COMMON.GEO'
1622 include 'COMMON.VAR'
1623 include 'COMMON.INTERACT'
1624 include 'COMMON.LOCAL'
1625 include 'COMMON.NAMES'
1626 include 'COMMON.CHAIN'
1627 include 'COMMON.FFIELD'
1628 include 'COMMON.SBRIDGE'
1629 include 'COMMON.HEADER'
1630 include 'COMMON.CONTROL'
1631 include 'COMMON.DBASE'
1632 include 'COMMON.THREAD'
1633 include 'COMMON.TIME1'
1635 C Set up variable list.
1641 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1643 ialph(i,1)=nvar+nside
1647 if (indphi.gt.0) then
1649 else if (indback.gt.0) then
1656 c----------------------------------------------------------------------------
1657 subroutine gen_dist_constr
1658 C Generate CA distance constraints.
1660 include 'DIMENSIONS'
1661 include 'COMMON.IOUNITS'
1662 include 'COMMON.GEO'
1663 include 'COMMON.VAR'
1664 include 'COMMON.INTERACT'
1665 include 'COMMON.LOCAL'
1666 include 'COMMON.NAMES'
1667 include 'COMMON.CHAIN'
1668 include 'COMMON.FFIELD'
1669 include 'COMMON.SBRIDGE'
1670 include 'COMMON.HEADER'
1671 include 'COMMON.CONTROL'
1672 include 'COMMON.DBASE'
1673 include 'COMMON.THREAD'
1674 include 'COMMON.SPLITELE'
1675 include 'COMMON.TIME1'
1676 integer i,j,itype_pdb(maxres)
1677 common /pizda/ itype_pdb
1679 double precision dist
1681 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1682 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1683 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1685 do i=nstart_sup,nstart_sup+nsup-1
1686 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1687 cd & ' seq_pdb', restyp(itype_pdb(i))
1688 do j=i+2,nstart_sup+nsup-1
1689 c 5/24/2020 Adam: Cutoff included to reduce array size
1691 if (dd.gt.r_cut_int) cycle
1693 ihpb(nhpb)=i+nstart_seq-nstart_sup
1694 jhpb(nhpb)=j+nstart_seq-nstart_sup
1699 cd write (iout,'(a)') 'Distance constraints:'
1704 cd if (ii.gt.nres) then
1709 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1710 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1711 cd & dhpb(i),forcon(i)
1715 c----------------------------------------------------------------------------
1718 include 'DIMENSIONS'
1719 include 'COMMON.MAP'
1720 include 'COMMON.IOUNITS'
1722 character*3 angid(4) /'THE','PHI','ALP','OME'/
1723 character*80 mapcard,ucase
1725 read (inp,'(a)') mapcard
1726 mapcard=ucase(mapcard)
1727 if (index(mapcard,'PHI').gt.0) then
1729 else if (index(mapcard,'THE').gt.0) then
1731 else if (index(mapcard,'ALP').gt.0) then
1733 else if (index(mapcard,'OME').gt.0) then
1736 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1737 stop 'Error - illegal variable spec in MAP card.'
1739 call readi (mapcard,'RES1',res1(imap),0)
1740 call readi (mapcard,'RES2',res2(imap),0)
1741 if (res1(imap).eq.0) then
1742 res1(imap)=res2(imap)
1743 else if (res2(imap).eq.0) then
1744 res2(imap)=res1(imap)
1746 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1748 & 'Error - illegal definition of variable group in MAP.'
1749 stop 'Error - illegal definition of variable group in MAP.'
1751 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1752 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1753 call readi(mapcard,'NSTEP',nstep(imap),0)
1754 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1756 & 'Illegal boundary and/or step size specification in MAP.'
1757 stop 'Illegal boundary and/or step size specification in MAP.'
1762 c----------------------------------------------------------------------------
1765 include 'DIMENSIONS'
1766 include 'COMMON.IOUNITS'
1767 include 'COMMON.GEO'
1768 include 'COMMON.CSA'
1769 include 'COMMON.BANK'
1770 include 'COMMON.CONTROL'
1772 character*620 mcmcard
1773 call card_concat(mcmcard)
1775 call readi(mcmcard,'NCONF',nconf,50)
1776 call readi(mcmcard,'NADD',nadd,0)
1777 call readi(mcmcard,'JSTART',jstart,1)
1778 call readi(mcmcard,'JEND',jend,1)
1779 call readi(mcmcard,'NSTMAX',nstmax,500000)
1780 call readi(mcmcard,'N0',n0,1)
1781 call readi(mcmcard,'N1',n1,6)
1782 call readi(mcmcard,'N2',n2,4)
1783 call readi(mcmcard,'N3',n3,0)
1784 call readi(mcmcard,'N4',n4,0)
1785 call readi(mcmcard,'N5',n5,0)
1786 call readi(mcmcard,'N6',n6,10)
1787 call readi(mcmcard,'N7',n7,0)
1788 call readi(mcmcard,'N8',n8,0)
1789 call readi(mcmcard,'N9',n9,0)
1790 call readi(mcmcard,'N14',n14,0)
1791 call readi(mcmcard,'N15',n15,0)
1792 call readi(mcmcard,'N16',n16,0)
1793 call readi(mcmcard,'N17',n17,0)
1794 call readi(mcmcard,'N18',n18,0)
1796 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1798 call readi(mcmcard,'NDIFF',ndiff,2)
1799 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1800 call readi(mcmcard,'IS1',is1,1)
1801 call readi(mcmcard,'IS2',is2,8)
1802 call readi(mcmcard,'NRAN0',nran0,4)
1803 call readi(mcmcard,'NRAN1',nran1,2)
1804 call readi(mcmcard,'IRR',irr,1)
1805 call readi(mcmcard,'NSEED',nseed,20)
1806 call readi(mcmcard,'NTOTAL',ntotal,10000)
1807 call reada(mcmcard,'CUT1',cut1,2.0d0)
1808 call reada(mcmcard,'CUT2',cut2,5.0d0)
1809 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1810 call readi(mcmcard,'ICMAX',icmax,3)
1811 call readi(mcmcard,'IRESTART',irestart,0)
1812 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1815 call reada(mcmcard,'DELE',dele,20.0d0)
1816 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1817 call readi(mcmcard,'IREF',iref,0)
1818 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1819 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1820 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1821 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1822 write (iout,*) "NCONF_IN",nconf_in
1825 c----------------------------------------------------------------------------
1828 include 'DIMENSIONS'
1829 include 'COMMON.MCM'
1830 include 'COMMON.MCE'
1831 include 'COMMON.IOUNITS'
1833 character*320 mcmcard
1835 call card_concat(mcmcard)
1836 call readi(mcmcard,'MAXACC',maxacc,100)
1837 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1838 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1839 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1840 call readi(mcmcard,'MAXREPM',maxrepm,200)
1841 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1842 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1843 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1844 call reada(mcmcard,'E_UP',e_up,5.0D0)
1845 call reada(mcmcard,'DELTE',delte,0.1D0)
1846 call readi(mcmcard,'NSWEEP',nsweep,5)
1847 call readi(mcmcard,'NSTEPH',nsteph,0)
1848 call readi(mcmcard,'NSTEPC',nstepc,0)
1849 call reada(mcmcard,'TMIN',tmin,298.0D0)
1850 call reada(mcmcard,'TMAX',tmax,298.0D0)
1851 call readi(mcmcard,'NWINDOW',nwindow,0)
1852 call readi(mcmcard,'PRINT_MC',print_mc,0)
1853 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1854 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1855 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1856 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1857 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1858 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1859 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1860 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1861 if (nwindow.gt.0) then
1862 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1864 winlen(i)=winend(i)-winstart(i)+1
1867 if (tmax.lt.tmin) tmax=tmin
1868 if (tmax.eq.tmin) then
1872 if (nstepc.gt.0 .and. nsteph.gt.0) then
1873 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1874 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1876 C Probabilities of different move types
1877 sumpro_type(0)=0.0D0
1878 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1879 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1880 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1881 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1882 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1883 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1884 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1886 print *,'i',i,' sumprotype',sumpro_type(i)
1887 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1888 print *,'i',i,' sumprotype',sumpro_type(i)
1892 c----------------------------------------------------------------------------
1893 subroutine read_minim
1895 include 'DIMENSIONS'
1896 include 'COMMON.MINIM'
1897 include 'COMMON.IOUNITS'
1899 character*320 minimcard
1900 call card_concat(minimcard)
1901 call readi(minimcard,'MAXMIN',maxmin,2000)
1902 call readi(minimcard,'MAXFUN',maxfun,5000)
1903 call readi(minimcard,'MINMIN',minmin,maxmin)
1904 call readi(minimcard,'MINFUN',minfun,maxmin)
1905 call reada(minimcard,'TOLF',tolf,1.0D-2)
1906 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1907 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1908 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1909 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1910 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1911 & 'Options in energy minimization:'
1912 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1913 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1914 & 'MinMin:',MinMin,' MinFun:',MinFun,
1915 & ' TolF:',TolF,' RTolF:',RTolF
1918 c----------------------------------------------------------------------------
1919 subroutine read_angles(kanal,*)
1921 include 'DIMENSIONS'
1922 include 'COMMON.GEO'
1923 include 'COMMON.VAR'
1924 include 'COMMON.CHAIN'
1925 include 'COMMON.IOUNITS'
1926 include 'COMMON.CONTROL'
1928 c Read angles from input
1930 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1931 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1932 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1933 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1936 c 9/7/01 avoid 180 deg valence angle
1937 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1939 theta(i)=deg2rad*theta(i)
1940 phi(i)=deg2rad*phi(i)
1941 alph(i)=deg2rad*alph(i)
1942 omeg(i)=deg2rad*omeg(i)
1947 c----------------------------------------------------------------------------
1948 subroutine reada(rekord,lancuch,wartosc,default)
1950 character*(*) rekord,lancuch
1951 double precision wartosc,default
1954 iread=index(rekord,lancuch)
1955 if (iread.eq.0) then
1959 iread=iread+ilen(lancuch)+1
1960 read (rekord(iread:),*,err=10,end=10) wartosc
1965 c----------------------------------------------------------------------------
1966 subroutine readi(rekord,lancuch,wartosc,default)
1968 character*(*) rekord,lancuch
1969 integer wartosc,default
1972 iread=index(rekord,lancuch)
1973 if (iread.eq.0) then
1977 iread=iread+ilen(lancuch)+1
1978 read (rekord(iread:),*,err=10,end=10) wartosc
1983 c----------------------------------------------------------------------------
1984 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1987 integer tablica(dim),default
1988 character*(*) rekord,lancuch
1995 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1996 if (iread.eq.0) return
1997 iread=iread+ilen(lancuch)+1
1998 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2001 c----------------------------------------------------------------------------
2002 subroutine multreada(rekord,lancuch,tablica,dim,default)
2005 double precision tablica(dim),default
2006 character*(*) rekord,lancuch
2013 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2014 if (iread.eq.0) return
2015 iread=iread+ilen(lancuch)+1
2016 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2019 c----------------------------------------------------------------------------
2020 subroutine openunits
2022 include 'DIMENSIONS'
2026 character*16 form,nodename
2029 include 'COMMON.SETUP'
2030 include 'COMMON.IOUNITS'
2032 include 'COMMON.CONTROL'
2033 integer lenpre,lenpot,ilen,lentmp,npos
2035 character*3 out1file_text,ucase
2040 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2041 call getenv_loc("PREFIX",prefix)
2043 call getenv_loc("POT",pot)
2044 call getenv_loc("DIRTMP",tmpdir)
2045 call getenv_loc("CURDIR",curdir)
2046 call getenv_loc("OUT1FILE",out1file_text)
2047 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2048 out1file_text=ucase(out1file_text)
2049 if (out1file_text(1:1).eq."Y") then
2052 out1file=fg_rank.gt.0
2057 if (lentmp.gt.0) then
2058 write (*,'(80(1h!))')
2059 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2060 write (*,'(80(1h!))')
2061 write (*,*)"All output files will be on node /tmp directory."
2063 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2064 if (me.eq.king) then
2065 write (*,*) "The master node is ",nodename
2066 else if (fg_rank.eq.0) then
2067 write (*,*) "I am the CG slave node ",nodename
2069 write (*,*) "I am the FG slave node ",nodename
2072 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2073 lenpre = lentmp+lenpre+1
2075 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2076 C Get the names and open the input files
2077 #if defined(WINIFL) || defined(WINPGI)
2078 open(1,file=pref_orig(:ilen(pref_orig))//
2079 & '.inp',status='old',readonly,shared)
2080 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2081 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2082 C Get parameter filenames and open the parameter files.
2083 call getenv_loc('BONDPAR',bondname)
2084 open (ibond,file=bondname,status='old',readonly,shared)
2085 call getenv_loc('THETPAR',thetname)
2086 open (ithep,file=thetname,status='old',readonly,shared)
2087 call getenv_loc('ROTPAR',rotname)
2088 open (irotam,file=rotname,status='old',readonly,shared)
2089 call getenv_loc('TORPAR',torname)
2090 open (itorp,file=torname,status='old',readonly,shared)
2091 call getenv_loc('TORDPAR',tordname)
2092 open (itordp,file=tordname,status='old',readonly,shared)
2093 call getenv_loc('FOURIER',fouriername)
2094 open (ifourier,file=fouriername,status='old',readonly,shared)
2095 call getenv_loc('ELEPAR',elename)
2096 open (ielep,file=elename,status='old',readonly,shared)
2097 call getenv_loc('SIDEPAR',sidename)
2098 open (isidep,file=sidename,status='old',readonly,shared)
2099 call getenv_loc('LIPTRANPAR',liptranname)
2100 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2101 #elif (defined CRAY) || (defined AIX)
2102 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2104 c print *,"Processor",myrank," opened file 1"
2105 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2106 c print *,"Processor",myrank," opened file 9"
2107 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2108 C Get parameter filenames and open the parameter files.
2109 call getenv_loc('BONDPAR',bondname)
2110 open (ibond,file=bondname,status='old',action='read')
2111 c print *,"Processor",myrank," opened file IBOND"
2112 call getenv_loc('THETPAR',thetname)
2113 open (ithep,file=thetname,status='old',action='read')
2115 call getenv_loc('THETPARPDB',thetname_pdb)
2116 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2118 c print *,"Processor",myrank," opened file ITHEP"
2119 call getenv_loc('ROTPAR',rotname)
2120 open (irotam,file=rotname,status='old',action='read')
2122 call getenv_loc('ROTPARPDB',rotname_pdb)
2123 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2125 c print *,"Processor",myrank," opened file IROTAM"
2126 call getenv_loc('TORPAR',torname)
2127 open (itorp,file=torname,status='old',action='read')
2128 c print *,"Processor",myrank," opened file ITORP"
2129 call getenv_loc('TORDPAR',tordname)
2130 open (itordp,file=tordname,status='old',action='read')
2131 c print *,"Processor",myrank," opened file ITORDP"
2132 call getenv_loc('SCCORPAR',sccorname)
2133 open (isccor,file=sccorname,status='old',action='read')
2134 c print *,"Processor",myrank," opened file ISCCOR"
2135 call getenv_loc('FOURIER',fouriername)
2136 open (ifourier,file=fouriername,status='old',action='read')
2137 c print *,"Processor",myrank," opened file IFOURIER"
2138 call getenv_loc('ELEPAR',elename)
2139 open (ielep,file=elename,status='old',action='read')
2140 c print *,"Processor",myrank," opened file IELEP"
2141 call getenv_loc('SIDEPAR',sidename)
2142 open (isidep,file=sidename,status='old',action='read')
2143 call getenv_loc('LIPTRANPAR',liptranname)
2144 open (iliptranpar,file=liptranname,status='old',action='read')
2145 c print *,"Processor",myrank," opened file ISIDEP"
2146 c print *,"Processor",myrank," opened parameter files"
2148 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2149 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2150 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2151 C Get parameter filenames and open the parameter files.
2152 call getenv_loc('BONDPAR',bondname)
2153 open (ibond,file=bondname,status='old')
2154 call getenv_loc('THETPAR',thetname)
2155 open (ithep,file=thetname,status='old')
2157 call getenv_loc('THETPARPDB',thetname_pdb)
2158 open (ithep_pdb,file=thetname_pdb,status='old')
2160 call getenv_loc('ROTPAR',rotname)
2161 open (irotam,file=rotname,status='old')
2163 call getenv_loc('ROTPARPDB',rotname_pdb)
2164 open (irotam_pdb,file=rotname_pdb,status='old')
2166 call getenv_loc('TORPAR',torname)
2167 open (itorp,file=torname,status='old')
2168 call getenv_loc('TORDPAR',tordname)
2169 open (itordp,file=tordname,status='old')
2170 call getenv_loc('SCCORPAR',sccorname)
2171 open (isccor,file=sccorname,status='old')
2172 call getenv_loc('FOURIER',fouriername)
2173 open (ifourier,file=fouriername,status='old')
2174 call getenv_loc('ELEPAR',elename)
2175 open (ielep,file=elename,status='old')
2176 call getenv_loc('SIDEPAR',sidename)
2177 open (isidep,file=sidename,status='old')
2178 call getenv_loc('LIPTRANPAR',liptranname)
2179 open (iliptranpar,file=liptranname,status='old')
2181 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2183 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2184 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2185 C Get parameter filenames and open the parameter files.
2186 call getenv_loc('BONDPAR',bondname)
2187 open (ibond,file=bondname,status='old',readonly)
2188 call getenv_loc('THETPAR',thetname)
2189 open (ithep,file=thetname,status='old',readonly)
2190 call getenv_loc('ROTPAR',rotname)
2191 open (irotam,file=rotname,status='old',readonly)
2192 call getenv_loc('TORPAR',torname)
2193 open (itorp,file=torname,status='old',readonly)
2194 call getenv_loc('TORDPAR',tordname)
2195 open (itordp,file=tordname,status='old',readonly)
2196 call getenv_loc('SCCORPAR',sccorname)
2197 open (isccor,file=sccorname,status='old',readonly)
2199 call getenv_loc('THETPARPDB',thetname_pdb)
2200 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2202 call getenv_loc('FOURIER',fouriername)
2203 open (ifourier,file=fouriername,status='old',readonly)
2204 call getenv_loc('ELEPAR',elename)
2205 open (ielep,file=elename,status='old',readonly)
2206 call getenv_loc('SIDEPAR',sidename)
2207 open (isidep,file=sidename,status='old',readonly)
2208 call getenv_loc('LIPTRANPAR',liptranname)
2209 open (iliptranpar,file=liptranname,status='old',action='read')
2211 call getenv_loc('ROTPARPDB',rotname_pdb)
2212 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2216 call getenv_loc('TUBEPAR',tubename)
2217 #if defined(WINIFL) || defined(WINPGI)
2218 open (itube,file=tubename,status='old',readonly,shared)
2219 #elif (defined CRAY) || (defined AIX)
2220 open (itube,file=tubename,status='old',action='read')
2222 open (itube,file=tubename,status='old')
2224 open (itube,file=tubename,status='old',readonly)
2229 C 8/9/01 In the newest version SCp interaction constants are read from a file
2230 C Use -DOLDSCP to use hard-coded constants instead.
2232 call getenv_loc('SCPPAR',scpname)
2233 #if defined(WINIFL) || defined(WINPGI)
2234 open (iscpp,file=scpname,status='old',readonly,shared)
2235 #elif (defined CRAY) || (defined AIX)
2236 open (iscpp,file=scpname,status='old',action='read')
2238 open (iscpp,file=scpname,status='old')
2240 open (iscpp,file=scpname,status='old',readonly)
2243 call getenv_loc('PATTERN',patname)
2244 #if defined(WINIFL) || defined(WINPGI)
2245 open (icbase,file=patname,status='old',readonly,shared)
2246 #elif (defined CRAY) || (defined AIX)
2247 open (icbase,file=patname,status='old',action='read')
2249 open (icbase,file=patname,status='old')
2251 open (icbase,file=patname,status='old',readonly)
2254 C Open output file only for CG processes
2255 c print *,"Processor",myrank," fg_rank",fg_rank
2256 if (fg_rank.eq.0) then
2258 if (nodes.eq.1) then
2261 npos = dlog10(dfloat(nodes-1))+1
2263 if (npos.lt.3) npos=3
2264 write (liczba,'(i1)') npos
2265 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2267 write (liczba,form) me
2268 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2269 & liczba(:ilen(liczba))
2270 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2272 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2274 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2275 & liczba(:ilen(liczba))//'.mol2'
2276 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2277 & liczba(:ilen(liczba))//'.stat'
2279 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2280 & //liczba(:ilen(liczba))//'.stat')
2281 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2284 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2285 & liczba(:ilen(liczba))//'.const'
2290 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2291 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2292 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2293 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2294 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2296 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2298 rest2name=prefix(:ilen(prefix))//'.rst'
2300 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2303 #if defined(AIX) || defined(PGI) || defined(CRAY)
2304 if (me.eq.king .or. .not. out1file) then
2305 open(iout,file=outname,status='unknown')
2307 open(iout,file="/dev/null",status="unknown")
2311 if (fg_rank.gt.0) then
2312 write (liczba,'(i3.3)') myrank/nfgtasks
2313 write (ll,'(bz,i3.3)') fg_rank
2314 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2320 open(igeom,file=intname,status='unknown',position='append')
2321 open(ipdb,file=pdbname,status='unknown')
2322 open(imol2,file=mol2name,status='unknown')
2323 open(istat,file=statname,status='unknown',position='append')
2325 c1out open(iout,file=outname,status='unknown')
2328 if (me.eq.king .or. .not.out1file)
2329 & open(iout,file=outname,status='unknown')
2331 if (fg_rank.gt.0) then
2332 write (liczba,'(i3.3)') myrank/nfgtasks
2333 write (ll,'(bz,i3.3)') fg_rank
2334 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2339 open(igeom,file=intname,status='unknown',access='append')
2340 open(ipdb,file=pdbname,status='unknown')
2341 open(imol2,file=mol2name,status='unknown')
2342 open(istat,file=statname,status='unknown',access='append')
2344 c1out open(iout,file=outname,status='unknown')
2347 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2348 csa_seed=prefix(:lenpre)//'.CSA.seed'
2349 csa_history=prefix(:lenpre)//'.CSA.history'
2350 csa_bank=prefix(:lenpre)//'.CSA.bank'
2351 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2352 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2353 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2354 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2355 csa_int=prefix(:lenpre)//'.int'
2356 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2357 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2358 csa_in=prefix(:lenpre)//'.CSA.in'
2359 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2362 write (iout,'(80(1h-))')
2363 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2364 write (iout,'(80(1h-))')
2365 write (iout,*) "Input file : ",
2366 & pref_orig(:ilen(pref_orig))//'.inp'
2367 write (iout,*) "Output file : ",
2368 & outname(:ilen(outname))
2370 write (iout,*) "Sidechain potential file : ",
2371 & sidename(:ilen(sidename))
2373 write (iout,*) "SCp potential file : ",
2374 & scpname(:ilen(scpname))
2376 write (iout,*) "Electrostatic potential file : ",
2377 & elename(:ilen(elename))
2378 write (iout,*) "Cumulant coefficient file : ",
2379 & fouriername(:ilen(fouriername))
2380 write (iout,*) "Torsional parameter file : ",
2381 & torname(:ilen(torname))
2382 write (iout,*) "Double torsional parameter file : ",
2383 & tordname(:ilen(tordname))
2384 write (iout,*) "SCCOR parameter file : ",
2385 & sccorname(:ilen(sccorname))
2386 write (iout,*) "Bond & inertia constant file : ",
2387 & bondname(:ilen(bondname))
2388 write (iout,*) "Bending-torsion parameter file : ",
2389 & thetname(:ilen(thetname))
2390 write (iout,*) "Rotamer parameter file : ",
2391 & rotname(:ilen(rotname))
2392 write (iout,*) "Thetpdb parameter file : ",
2393 & thetname_pdb(:ilen(thetname_pdb))
2394 write (iout,*) "Threading database : ",
2395 & patname(:ilen(patname))
2397 &write (iout,*)" DIRTMP : ",
2399 write (iout,'(80(1h-))')
2403 c----------------------------------------------------------------------------
2404 subroutine card_concat(card)
2405 implicit real*8 (a-h,o-z)
2406 include 'DIMENSIONS'
2407 include 'COMMON.IOUNITS'
2409 character*80 karta,ucase
2411 read (inp,'(a)') karta
2414 do while (karta(80:80).eq.'&')
2415 card=card(:ilen(card)+1)//karta(:79)
2416 read (inp,'(a)') karta
2419 card=card(:ilen(card)+1)//karta
2422 c------------------------------------------------------------------------------
2425 include 'DIMENSIONS'
2426 include 'COMMON.CHAIN'
2427 include 'COMMON.IOUNITS'
2428 include 'COMMON.CONTROL'
2430 include 'COMMON.QRESTR'
2432 include 'COMMON.LAGRANGE.5diag'
2434 include 'COMMON.LAGRANGE'
2437 open(irest2,file=rest2name,status='unknown')
2438 read(irest2,*) totT,EK,potE,totE,t_bath
2441 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2444 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2447 read (irest2,*) iset
2452 c------------------------------------------------------------------------------
2453 subroutine read_fragments
2455 include 'DIMENSIONS'
2459 include 'COMMON.SETUP'
2460 include 'COMMON.CHAIN'
2461 include 'COMMON.IOUNITS'
2463 include 'COMMON.QRESTR'
2464 include 'COMMON.CONTROL'
2466 read(inp,*) nset,nfrag,npair,nfrag_back
2467 loc_qlike=(nfrag_back.lt.0)
2468 nfrag_back=iabs(nfrag_back)
2469 if(me.eq.king.or..not.out1file)
2470 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2471 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2473 read(inp,*) mset(iset)
2475 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2477 if(me.eq.king.or..not.out1file)
2478 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2479 & ifrag(2,i,iset), qinfrag(i,iset)
2482 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2484 if(me.eq.king.or..not.out1file)
2485 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2486 & ipair(2,i,iset), qinpair(i,iset)
2490 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2491 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2492 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2493 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2494 if(me.eq.king.or..not.out1file)
2495 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2496 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2497 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2498 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2502 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2503 & wfrag_back(3,i,iset),
2504 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2505 if(me.eq.king.or..not.out1file)
2506 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2507 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2513 C---------------------------------------------------------------------------
2514 subroutine read_afminp
2516 include 'DIMENSIONS'
2520 include 'COMMON.SETUP'
2521 include 'COMMON.CONTROL'
2522 include 'COMMON.CHAIN'
2523 include 'COMMON.IOUNITS'
2524 include 'COMMON.SBRIDGE'
2525 character*320 afmcard
2527 c print *, "wchodze"
2528 call card_concat(afmcard)
2529 call readi(afmcard,"BEG",afmbeg,0)
2530 call readi(afmcard,"END",afmend,0)
2531 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2532 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2533 c print *,'FORCE=' ,forceAFMconst
2534 CCCC NOW PROPERTIES FOR AFM
2537 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2539 distafminit=dsqrt(distafminit)
2540 c print *,'initdist',distafminit
2543 c-------------------------------------------------------------------------------
2544 subroutine read_saxs_constr
2546 include 'DIMENSIONS'
2550 include 'COMMON.SETUP'
2551 include 'COMMON.CONTROL'
2552 include 'COMMON.SAXS'
2553 include 'COMMON.CHAIN'
2554 include 'COMMON.IOUNITS'
2555 include 'COMMON.SBRIDGE'
2556 double precision cm(3),cnorm
2559 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2561 if (saxs_mode.eq.0) then
2562 c SAXS distance distribution
2564 read(inp,*) distsaxs(i),Psaxs(i)
2568 Cnorm = Cnorm + Psaxs(i)
2570 write (iout,*) "Cnorm",Cnorm
2572 Psaxs(i)=Psaxs(i)/Cnorm
2574 write (iout,*) "Normalized distance distribution from SAXS"
2576 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2580 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2582 write (iout,*) "Wsaxs0",Wsaxs0
2586 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2593 cm(j)=cm(j)+Csaxs(j,i)
2601 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2604 write (iout,*) "SAXS sphere coordinates"
2606 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2612 c-------------------------------------------------------------------------------
2613 subroutine read_dist_constr
2615 include 'DIMENSIONS'
2619 include 'COMMON.SETUP'
2620 include 'COMMON.CONTROL'
2621 include 'COMMON.CHAIN'
2622 include 'COMMON.IOUNITS'
2623 include 'COMMON.SBRIDGE'
2624 include 'COMMON.INTERACT'
2625 integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2626 integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2627 double precision wfrag_(100),wpair_(1000)
2628 double precision ddjk,dist,dist_cut,fordepthmax
2629 character*5000 controlcard
2630 logical normalize,next
2632 double precision scal_bfac
2633 double precision xlink(4,0:4) /
2635 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2636 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2637 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2638 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2639 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2640 c print *, "WCHODZE"
2641 write (iout,*) "Calling read_dist_constr"
2642 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2644 restr_on_coord=.false.
2653 call card_concat(controlcard)
2654 next = index(controlcard,"NEXT").gt.0
2655 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2656 write (iout,*) "restr_type",restr_type
2657 call readi(controlcard,"NFRAG",nfrag_,0)
2658 call readi(controlcard,"NFRAG",nfrag_,0)
2659 call readi(controlcard,"NPAIR",npair_,0)
2660 call readi(controlcard,"NDIST",ndist_,0)
2661 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2662 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2663 if (restr_type.eq.10)
2664 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2665 if (restr_type.eq.12)
2666 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2667 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2668 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2669 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2670 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2671 normalize = index(controlcard,"NORMALIZE").gt.0
2672 write (iout,*) "WBOLTZD",wboltzd
2673 write (iout,*) "SCAL_PEAK",scal_peak
2674 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2675 write (iout,*) "IFRAG"
2677 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2679 write (iout,*) "IPAIR"
2681 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2683 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2685 & "Distance restraints as generated from reference structure"
2687 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2688 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2689 & ifrag_(2,i)=nstart_sup+nsup-1
2690 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2692 if (wfrag_(i).eq.0.0d0) cycle
2693 do j=ifrag_(1,i),ifrag_(2,i)-1
2694 do k=j+1,ifrag_(2,i)
2695 c write (iout,*) "j",j," k",k
2697 if (restr_type.eq.1) then
2703 forcon(nhpb)=wfrag_(i)
2704 else if (constr_dist.eq.2) then
2705 if (ddjk.le.dist_cut) then
2711 forcon(nhpb)=wfrag_(i)
2713 else if (restr_type.eq.3) then
2719 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2722 if (.not.out1file .or. me.eq.king)
2723 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2724 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2726 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2727 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2733 if (wpair_(i).eq.0.0d0) cycle
2741 do j=ifrag_(1,ii),ifrag_(2,ii)
2742 do k=ifrag_(1,jj),ifrag_(2,jj)
2744 if (restr_type.eq.1) then
2750 forcon(nhpb)=wpair_(i)
2751 else if (constr_dist.eq.2) then
2752 if (ddjk.le.dist_cut) then
2758 forcon(nhpb)=wpair_(i)
2760 else if (restr_type.eq.3) then
2766 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2769 if (.not.out1file .or. me.eq.king)
2770 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2771 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2773 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2774 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2781 write (iout,*) "Distance restraints as read from input"
2783 if (restr_type.eq.12) then
2784 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2785 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2786 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2787 & fordepth_peak(nhpb_peak+1),npeak
2788 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2789 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2790 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2791 c & fordepth_peak(nhpb_peak+1),npeak
2792 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2793 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2794 nhpb_peak=nhpb_peak+1
2795 irestr_type_peak(nhpb_peak)=12
2796 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2799 if (.not.out1file .or. me.eq.king)
2800 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2801 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2802 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2803 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2804 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2806 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2807 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2808 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2809 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2810 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2812 if (ibecarb_peak(nhpb_peak).eq.3) then
2813 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2814 else if (ibecarb_peak(nhpb_peak).eq.2) then
2815 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2816 else if (ibecarb_peak(nhpb_peak).eq.1) then
2817 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2818 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2820 else if (restr_type.eq.11) then
2821 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2822 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2823 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2824 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2826 irestr_type(nhpb)=11
2828 if (.not.out1file .or. me.eq.king)
2829 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2830 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2831 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2833 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2834 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2835 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2837 c if (ibecarb(nhpb).gt.0) then
2838 c ihpb(nhpb)=ihpb(nhpb)+nres
2839 c jhpb(nhpb)=jhpb(nhpb)+nres
2841 if (ibecarb(nhpb).eq.3) then
2842 ihpb(nhpb)=ihpb(nhpb)+nres
2843 else if (ibecarb(nhpb).eq.2) then
2844 ihpb(nhpb)=ihpb(nhpb)+nres
2845 else if (ibecarb(nhpb).eq.1) then
2846 ihpb(nhpb)=ihpb(nhpb)+nres
2847 jhpb(nhpb)=jhpb(nhpb)+nres
2849 else if (restr_type.eq.10) then
2850 c Cross-lonk Markov-like potential
2851 call card_concat(controlcard)
2852 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2853 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2855 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2856 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2857 if (index(controlcard,"ZL").gt.0) then
2859 else if (index(controlcard,"ADH").gt.0) then
2861 else if (index(controlcard,"PDH").gt.0) then
2863 else if (index(controlcard,"DSS").gt.0) then
2868 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2869 & xlink(1,link_type))
2870 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2871 & xlink(2,link_type))
2872 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2873 & xlink(3,link_type))
2874 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2875 & xlink(4,link_type))
2876 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2877 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2878 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2879 if (forcon(nhpb+1).le.0.0d0 .or.
2880 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2882 irestr_type(nhpb)=10
2883 if (ibecarb(nhpb).eq.3) then
2884 jhpb(nhpb)=jhpb(nhpb)+nres
2885 else if (ibecarb(nhpb).eq.2) then
2886 ihpb(nhpb)=ihpb(nhpb)+nres
2887 else if (ibecarb(nhpb).eq.1) then
2888 ihpb(nhpb)=ihpb(nhpb)+nres
2889 jhpb(nhpb)=jhpb(nhpb)+nres
2892 if (.not.out1file .or. me.eq.king)
2893 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2894 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2895 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2898 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2899 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2900 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2905 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2906 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2907 if (forcon(nhpb+1).gt.0.0d0) then
2909 if (dhpb1(nhpb).eq.0.0d0) then
2914 if (ibecarb(nhpb).eq.3) then
2915 jhpb(nhpb)=jhpb(nhpb)+nres
2916 else if (ibecarb(nhpb).eq.2) then
2917 ihpb(nhpb)=ihpb(nhpb)+nres
2918 else if (ibecarb(nhpb).eq.1) then
2919 ihpb(nhpb)=ihpb(nhpb)+nres
2920 jhpb(nhpb)=jhpb(nhpb)+nres
2922 if (dhpb(nhpb).eq.0.0d0)
2923 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2926 if (.not.out1file .or. me.eq.king)
2927 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2928 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2930 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2931 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2934 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2935 C if (forcon(nhpb+1).gt.0.0d0) then
2937 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2940 if (restr_type.eq.4) then
2941 write (iout,*) "The BFAC array"
2943 write (iout,'(i5,f10.5)') i,bfac(i)
2946 if (itype(i).eq.ntyp1) cycle
2948 if (itype(j).eq.ntyp1) cycle
2949 if (itype(i).eq.10) then
2954 if (itype(j).eq.10) then
2964 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2967 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2968 ihpb(nhpb)=i+nres*ii
2969 jhpb(nhpb)=j+nres*jj
2970 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2972 if (.not.out1file .or. me.eq.king) then
2973 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2974 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2975 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2979 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2980 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2981 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2991 if (restr_type.eq.5) then
2992 restr_on_coord=.true.
2994 if (itype(i).eq.ntyp1) cycle
2995 bfac(i)=(scal_bfac/bfac(i))**2
3004 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3005 & fordepthmax=fordepth(i)
3008 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3013 c-------------------------------------------------------------------------------
3014 subroutine read_constr_homology
3016 include 'DIMENSIONS'
3020 include 'COMMON.SETUP'
3021 include 'COMMON.CONTROL'
3022 include 'COMMON.HOMOLOGY'
3023 include 'COMMON.CHAIN'
3024 include 'COMMON.IOUNITS'
3026 include 'COMMON.QRESTR'
3027 include 'COMMON.GEO'
3028 include 'COMMON.INTERACT'
3029 include 'COMMON.NAMES'
3031 c For new homol impl
3033 include 'COMMON.VAR'
3036 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3038 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3039 c & sigma_odl_temp(maxres,maxres,max_template)
3041 character*24 model_ki_dist, model_ki_angle
3042 character*500 controlcard
3043 integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
3044 & ik,iistart,nres_temp
3047 logical liiflag,lfirst
3050 c FP - Nov. 2014 Temporary specifications for new vars
3052 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3054 double precision, dimension (max_template,maxres) :: rescore
3055 double precision, dimension (max_template,maxres) :: rescore2
3056 double precision, dimension (max_template,maxres) :: rescore3
3057 double precision distal
3058 character*24 tpl_k_rescore
3059 character*256 pdbfile
3060 c -----------------------------------------------------------------
3061 c Reading multiple PDB ref structures and calculation of retraints
3062 c not using pre-computed ones stored in files model_ki_{dist,angle}
3064 c -----------------------------------------------------------------
3067 c Alternative: reading from input
3068 call card_concat(controlcard)
3069 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3070 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3071 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3072 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3073 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3074 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3075 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3076 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3077 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3078 if(.not.read2sigma.and.start_from_model) then
3079 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3080 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3081 start_from_model=.false.
3082 iranconf=(indpdb.le.0)
3084 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3085 & write(iout,*) 'START_FROM_MODELS is ON'
3086 c if(start_from_model .and. rest) then
3087 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3088 c write(iout,*) 'START_FROM_MODELS is OFF'
3089 c write(iout,*) 'remove restart keyword from input'
3092 if (start_from_model) nmodel_start=constr_homology
3093 if (homol_nset.gt.1)then
3094 call card_concat(controlcard)
3095 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3096 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3097 write(iout,*) "iset homology_weight "
3099 write(iout,*) i,waga_homology(i)
3102 iset=mod(kolor,homol_nset)+1
3105 waga_homology(1)=1.0
3108 cd write (iout,*) "nnt",nnt," nct",nct
3115 c write(iout,*) 'nnt=',nnt,'nct=',nct
3118 do k=1,constr_homology
3131 if (read_homol_frag) then
3132 call read_klapaucjusz
3135 do k=1,constr_homology
3137 read(inp,'(a)') pdbfile
3138 pdbfiles_chomo(k)=pdbfile
3139 if(me.eq.king .or. .not. out1file)
3140 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3141 & pdbfile(:ilen(pdbfile))
3142 open(ipdbin,file=pdbfile,status='old',err=33)
3144 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3145 & pdbfile(:ilen(pdbfile))
3148 c print *,'Begin reading pdb data'
3150 c Files containing res sim or local scores (former containing sigmas)
3153 write(kic2,'(bz,i2.2)') k
3155 tpl_k_rescore="template"//kic2//".sco"
3159 if (read2sigma) then
3160 call readpdb_template(k)
3167 c Distance restraints
3170 C Copy the coordinates from reference coordinates (?)
3171 do i=1,2*nres_chomo(k)
3174 c write (iout,*) "c(",j,i,") =",c(j,i)
3178 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3180 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3181 open (ientin,file=tpl_k_rescore,status='old')
3182 if (nnt.gt.1) rescore(k,1)=0.0d0
3183 do irec=nnt,nct ! loop for reading res sim
3184 if (read2sigma) then
3185 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3186 & rescore3_tmp,idomain_tmp
3188 idomain(k,i_tmp)=idomain_tmp
3189 rescore(k,i_tmp)=rescore_tmp
3190 rescore2(k,i_tmp)=rescore2_tmp
3191 rescore3(k,i_tmp)=rescore3_tmp
3192 if (.not. out1file .or. me.eq.king)
3193 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3194 & i_tmp,rescore2_tmp,rescore_tmp,
3195 & rescore3_tmp,idomain_tmp
3198 read (ientin,*,end=1401) rescore_tmp
3200 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3201 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3202 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3207 if (waga_dist.ne.0.0d0) then
3215 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3216 c write (iout,*) k,i,j,distal,dist2_cut
3218 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3219 & .and. distal.le.dist2_cut ) then
3225 c write (iout,*) "k",k
3226 c write (iout,*) "i",i," j",j," constr_homology",
3231 if (read2sigma) then
3234 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3236 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3237 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3238 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3240 if (odl(k,ii).le.dist_cut) then
3241 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3244 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3245 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3247 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3248 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3252 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3255 l_homo(k,ii)=.false.
3261 c write (iout,*) "Distance restraints set"
3264 c Theta, dihedral and SC retraints
3266 if (waga_angle.gt.0.0d0) then
3267 c open (ientin,file=tpl_k_sigma_dih,status='old')
3268 c do irec=1,maxres-3 ! loop for reading sigma_dih
3269 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3270 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3271 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3272 c & sigma_dih(k,i+nnt-1)
3277 if (idomain(k,i).eq.0) then
3281 dih(k,i)=phiref(i) ! right?
3282 c read (ientin,*) sigma_dih(k,i) ! original variant
3283 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3284 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3285 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3286 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3287 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3289 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3290 & rescore(k,i-2)+rescore(k,i-3))/4.0
3291 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3292 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3293 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3294 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3295 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3296 c Instead of res sim other local measure of b/b str reliability possible
3297 if (sigma_dih(k,i).ne.0)
3298 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3299 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3303 c write (iout,*) "Dihedral angle restraints set"
3306 if (waga_theta.gt.0.0d0) then
3307 c open (ientin,file=tpl_k_sigma_theta,status='old')
3308 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3309 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3310 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3311 c & sigma_theta(k,i+nnt-1)
3316 do i = nnt+2,nct ! right? without parallel.
3317 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3318 c do i=ithet_start,ithet_end ! with FG parallel.
3319 if (idomain(k,i).eq.0) then
3320 sigma_theta(k,i)=0.0
3323 thetatpl(k,i)=thetaref(i)
3324 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3325 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3326 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3327 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3328 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3329 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3330 & rescore(k,i-2))/3.0
3331 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3332 if (sigma_theta(k,i).ne.0)
3333 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3335 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3336 c rescore(k,i-2) ! right expression ?
3337 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3340 c write (iout,*) "Angle restraints set"
3343 if (waga_d.gt.0.0d0) then
3344 c open (ientin,file=tpl_k_sigma_d,status='old')
3345 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3346 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3347 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3348 c & sigma_d(k,i+nnt-1)
3352 do i = nnt,nct ! right? without parallel.
3353 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3354 c do i=loc_start,loc_end ! with FG parallel.
3355 if (itype(i).eq.10) cycle
3356 if (idomain(k,i).eq.0 ) then
3363 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3364 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3365 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3366 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3367 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3368 if (sigma_d(k,i).ne.0)
3369 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3371 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3372 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3373 c read (ientin,*) sigma_d(k,i) ! 1st variant
3377 c write (iout,*) "SC restraints set"
3380 c remove distance restraints not used in any model from the list
3381 c shift data in all arrays
3383 c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
3384 if (waga_dist.ne.0.0d0) then
3391 c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3392 c & .and. distal.le.dist2_cut ) then
3393 c write (iout,*) "i",i," j",j," ii",ii
3395 if (ii_in_use(ii).eq.0.and.liiflag.or.
3396 & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3403 if(i10.eq.lim_odl) i10=i10+1
3405 ires_homo(iistart+ki)=ires_homo(ki+i01)
3406 jres_homo(iistart+ki)=jres_homo(ki+i01)
3407 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3408 do k=1,constr_homology
3409 odl(k,iistart+ki)=odl(k,ki+i01)
3410 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3411 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3414 iistart=iistart+i10-i01
3417 if (ii_in_use(ii).ne.0.and..not.liiflag) then
3425 c write (iout,*) "Removing distances completed"
3427 endif ! .not. klapaucjusz
3429 if (constr_homology.gt.0) call homology_partition
3430 c write (iout,*) "After homology_partition"
3432 if (constr_homology.gt.0) call init_int_table
3433 c write (iout,*) "After init_int_table"
3435 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3436 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3440 if (.not.out_template_restr) return
3441 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3442 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3443 write (iout,*) "Distance restraints from templates"
3445 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3446 & ii,ires_homo(ii),jres_homo(ii),
3447 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3448 & ki=1,constr_homology)
3450 write (iout,*) "Dihedral angle restraints from templates"
3452 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3453 & (rad2deg*dih(ki,i),
3454 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3456 write (iout,*) "Virtual-bond angle restraints from templates"
3458 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3459 & (rad2deg*thetatpl(ki,i),
3460 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3462 write (iout,*) "SC restraints from templates"
3464 write(iout,'(i5,100(4f8.2,4x))') i,
3465 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3466 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3469 c -----------------------------------------------------------------
3472 c----------------------------------------------------------------------
3474 subroutine flush(iu)
3479 subroutine flush(iu)
3484 c------------------------------------------------------------------------------
3485 subroutine copy_to_tmp(source)
3487 include "DIMENSIONS"
3488 include "COMMON.IOUNITS"
3489 character*(*) source
3490 character* 256 tmpfile
3494 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3495 inquire(file=tmpfile,exist=ex)
3497 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3498 & " to temporary directory..."
3499 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3500 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3504 c------------------------------------------------------------------------------
3505 subroutine move_from_tmp(source)
3507 include "DIMENSIONS"
3508 include "COMMON.IOUNITS"
3509 character*(*) source
3512 write (*,*) "Moving ",source(:ilen(source)),
3513 & " from temporary directory to working directory"
3514 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3515 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3518 c------------------------------------------------------------------------------
3519 subroutine random_init(seed)
3521 C Initialize random number generator
3524 include 'DIMENSIONS'
3527 logical OKRandom, prng_restart
3529 integer iseed_array(4)
3530 integer error_msg,ierr
3532 include 'COMMON.IOUNITS'
3533 include 'COMMON.TIME1'
3534 include 'COMMON.THREAD'
3535 include 'COMMON.SBRIDGE'
3536 include 'COMMON.CONTROL'
3537 include 'COMMON.MCM'
3538 include 'COMMON.MAP'
3539 include 'COMMON.HEADER'
3540 include 'COMMON.CSA'
3541 include 'COMMON.CHAIN'
3542 include 'COMMON.MUCA'
3544 include 'COMMON.FFIELD'
3545 include 'COMMON.SETUP'
3547 double precision seed,ran_number
3548 iseed=-dint(dabs(seed))
3549 if (iseed.eq.0) then
3550 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3551 & 'Random seed undefined. The program will stop.'
3552 write (*,'(/80(1h*)/20x,a/80(1h*))')
3553 & 'Random seed undefined. The program will stop.'
3555 call mpi_finalize(mpi_comm_world,ierr)
3557 stop 'Bad random seed.'
3560 if (fg_rank.eq.0) then
3564 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3565 OKRandom = prng_restart(me,iseed)
3568 tmp=65536.0d0**(4-i)
3569 iseed_array(i) = dint(seed/tmp)
3570 seed=seed-iseed_array(i)*tmp
3573 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3574 & (iseed_array(i),i=1,4)
3575 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3576 & (iseed_array(i),i=1,4)
3577 OKRandom = prng_restart(me,iseed_array)
3580 c r1 = prng_next(me)
3581 r1=ran_number(0.0D0,1.0D0)
3583 & write (iout,*) 'ran_num',r1
3584 if (r1.lt.0.0d0) OKRandom=.false.
3586 if (.not.OKRandom) then
3587 write (iout,*) 'PRNG IS NOT WORKING!!!'
3588 print *,'PRNG IS NOT WORKING!!!'
3591 call mpi_abort(mpi_comm_world,error_msg,ierr)
3594 write (iout,*) 'too many processors for parallel prng'
3595 write (*,*) 'too many processors for parallel prng'
3603 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3607 c----------------------------------------------------------------------
3608 subroutine read_klapaucjusz
3610 include 'DIMENSIONS'
3614 include 'COMMON.SETUP'
3615 include 'COMMON.CONTROL'
3616 include 'COMMON.HOMOLOGY'
3617 include 'COMMON.CHAIN'
3618 include 'COMMON.IOUNITS'
3620 include 'COMMON.GEO'
3621 include 'COMMON.INTERACT'
3622 include 'COMMON.NAMES'
3623 character*256 fragfile
3624 integer ninclust(maxclust),inclust(max_template,maxclust),
3625 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
3628 character*24 model_ki_dist, model_ki_angle
3629 character*500 controlcard
3630 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3631 & ik,ll,ii,kk,iistart,iishift,lim_xx
3632 double precision distal
3633 logical lprn /.true./
3640 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3641 double precision, dimension (max_template,maxres) :: rescore
3642 double precision, dimension (max_template,maxres) :: rescore2
3643 character*24 tpl_k_rescore
3644 character*256 pdbfile
3647 c For new homol impl
3649 include 'COMMON.VAR'
3651 call getenv("FRAGFILE",fragfile)
3652 open(ientin,file=fragfile,status="old",err=10)
3653 read(ientin,*) constr_homology,nclust
3659 do k=1,constr_homology
3660 read(ientin,'(a)') pdbfile
3661 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3662 & pdbfile(:ilen(pdbfile))
3663 pdbfiles_chomo(k)=pdbfile
3664 open(ipdbin,file=pdbfile,status='old',err=33)
3666 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3667 & pdbfile(:ilen(pdbfile))
3672 call readpdb_template(k)
3682 read(ientin,*) ninclust(i),nresclust(i)
3683 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3684 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3687 c Loop over clusters
3690 do ll = 1,ninclust(l)
3698 idomain(k,iresclust(i,l)+1) = 1
3700 idomain(k,iresclust(i,l)) = 1
3704 c Distance restraints
3707 C Copy the coordinates from reference coordinates (?)
3713 c write (iout,*) "c(",j,i,") =",c(j,i)
3716 call int_from_cart(.true.,.false.)
3717 call sc_loc_geom(.false.)
3719 thetaref(i)=theta(i)
3723 if (waga_dist.ne.0.0d0) then
3731 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3732 c write (iout,*) k,i,j,distal,dist2_cut
3734 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3735 & .and. distal.le.dist2_cut ) then
3741 c write (iout,*) "k",k
3742 c write (iout,*) "i",i," j",j," constr_homology",
3747 if (read2sigma) then
3750 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3752 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3753 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3754 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3756 if (odl(k,ii).le.dist_cut) then
3757 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3760 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3761 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3763 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3764 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3768 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3771 c l_homo(k,ii)=.false.
3778 c Theta, dihedral and SC retraints
3780 if (waga_angle.gt.0.0d0) then
3782 if (idomain(k,i).eq.0) then
3783 c sigma_dih(k,i)=0.0
3787 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3788 & rescore(k,i-2)+rescore(k,i-3))/4.0
3789 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3790 c & " sigma_dihed",sigma_dih(k,i)
3791 if (sigma_dih(k,i).ne.0)
3792 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3797 if (waga_theta.gt.0.0d0) then
3799 if (idomain(k,i).eq.0) then
3800 c sigma_theta(k,i)=0.0
3803 thetatpl(k,i)=thetaref(i)
3804 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3805 & rescore(k,i-2))/3.0
3806 if (sigma_theta(k,i).ne.0)
3807 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3811 if (waga_d.gt.0.0d0) then
3813 if (itype(i).eq.10) cycle
3814 if (idomain(k,i).eq.0 ) then
3821 sigma_d(k,i)=rescore(k,i)
3822 if (sigma_d(k,i).ne.0)
3823 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3824 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3830 c remove distance restraints not used in any model from the list
3831 c shift data in all arrays
3833 if (waga_dist.ne.0.0d0) then
3839 if (ii_in_use(ii).eq.0.and.liiflag) then
3843 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3844 & .not.liiflag.and.ii.eq.lim_odl) then
3845 if (ii.eq.lim_odl) then
3846 iishift=ii-iistart+1
3851 do ki=iistart,lim_odl-iishift
3852 ires_homo(ki)=ires_homo(ki+iishift)
3853 jres_homo(ki)=jres_homo(ki+iishift)
3854 ii_in_use(ki)=ii_in_use(ki+iishift)
3855 do k=1,constr_homology
3856 odl(k,ki)=odl(k,ki+iishift)
3857 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3858 l_homo(k,ki)=l_homo(k,ki+iishift)
3862 lim_odl=lim_odl-iishift
3869 10 stop "Error in fragment file"