2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
85 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
87 character*320 controlcard
92 read (INP,'(a)') titel
93 call card_concat(controlcard)
94 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
95 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
96 call reada(controlcard,'SEED',seed,0.0D0)
97 call random_init(seed)
98 C Set up the time limit (caution! The time must be input in minutes!)
99 read_cart=index(controlcard,'READ_CART').gt.0
100 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
101 write (iout,*) "constr_dist",constr_dist
102 call readi(controlcard,'NSAXS',nsaxs,0)
103 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
104 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
106 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
107 call readi(controlcard,'SYM',symetr,1)
108 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
109 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
110 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
111 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
112 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
113 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
114 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
115 call reada(controlcard,'DRMS',drms,0.1D0)
116 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
117 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
118 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
119 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
120 write (iout,'(a,f10.1)')'DRMS = ',drms
121 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
122 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
124 call readi(controlcard,'NZ_START',nz_start,0)
125 call readi(controlcard,'NZ_END',nz_end,0)
126 call readi(controlcard,'IZ_SC',iz_sc,0)
128 safety = 60.0d0*safety
131 call reada(controlcard,"T_BATH",t_bath,300.0d0)
132 minim=(index(controlcard,'MINIMIZE').gt.0)
133 dccart=(index(controlcard,'CART').gt.0)
134 overlapsc=(index(controlcard,'OVERLAP').gt.0)
135 overlapsc=.not.overlapsc
136 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
137 searchsc=.not.searchsc
138 sideadd=(index(controlcard,'SIDEADD').gt.0)
139 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
140 outpdb=(index(controlcard,'PDBOUT').gt.0)
141 outmol2=(index(controlcard,'MOL2OUT').gt.0)
142 pdbref=(index(controlcard,'PDBREF').gt.0)
143 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
144 indpdb=index(controlcard,'PDBSTART')
145 extconf=(index(controlcard,'EXTCONF').gt.0)
146 AFMlog=(index(controlcard,'AFM'))
147 selfguide=(index(controlcard,'SELFGUIDE'))
148 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
149 call readi(controlcard,'IPRINT',iprint,0)
150 call readi(controlcard,'MAXGEN',maxgen,10000)
151 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
152 call readi(controlcard,"KDIAG",kdiag,0)
153 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
154 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
155 & write (iout,*) "RESCALE_MODE",rescale_mode
156 split_ene=index(controlcard,'SPLIT_ENE').gt.0
157 if (index(controlcard,'REGULAR').gt.0.0D0) then
158 call reada(controlcard,'WEIDIS',weidis,0.1D0)
162 if (index(controlcard,'CHECKGRAD').gt.0) then
164 if (index(controlcard,'CART').gt.0) then
166 elseif (index(controlcard,'CARINT').gt.0) then
171 elseif (index(controlcard,'THREAD').gt.0) then
173 call readi(controlcard,'THREAD',nthread,0)
174 if (nthread.gt.0) then
175 call reada(controlcard,'WEIDIS',weidis,0.1D0)
178 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
179 stop 'Error termination in Read_Control.'
181 else if (index(controlcard,'MCMA').gt.0) then
183 else if (index(controlcard,'MCEE').gt.0) then
185 else if (index(controlcard,'MULTCONF').gt.0) then
187 else if (index(controlcard,'MAP').gt.0) then
189 call readi(controlcard,'MAP',nmap,0)
190 else if (index(controlcard,'CSA').gt.0) then
192 crc else if (index(controlcard,'ZSCORE').gt.0) then
194 crc ZSCORE is rm from UNRES, modecalc=9 is available
197 cfcm else if (index(controlcard,'MCMF').gt.0) then
199 else if (index(controlcard,'SOFTREG').gt.0) then
201 else if (index(controlcard,'CHECK_BOND').gt.0) then
203 else if (index(controlcard,'TEST').gt.0) then
205 else if (index(controlcard,'MD').gt.0) then
207 else if (index(controlcard,'RE ').gt.0) then
211 lmuca=index(controlcard,'MUCA').gt.0
212 call readi(controlcard,'MUCADYN',mucadyn,0)
213 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
214 if (lmuca .and. (me.eq.king .or. .not.out1file ))
216 write (iout,*) 'MUCADYN=',mucadyn
217 write (iout,*) 'MUCASMOOTH=',muca_smooth
220 iscode=index(controlcard,'ONE_LETTER')
221 indphi=index(controlcard,'PHI')
222 indback=index(controlcard,'BACK')
223 iranconf=index(controlcard,'RAND_CONF')
224 i2ndstr=index(controlcard,'USE_SEC_PRED')
225 gradout=index(controlcard,'GRADOUT').gt.0
226 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
227 C DISTCHAINMAX become obsolete for periodic boundry condition
228 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
229 C Reading the dimensions of box in x,y,z coordinates
230 call reada(controlcard,'BOXX',boxxsize,100.0d0)
231 call reada(controlcard,'BOXY',boxysize,100.0d0)
232 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
233 c Cutoff range for interactions
234 call reada(controlcard,"R_CUT",r_cut,15.0d0)
235 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
236 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
237 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
238 if (lipthick.gt.0.0d0) then
239 bordliptop=(boxzsize+lipthick)/2.0
240 bordlipbot=bordliptop-lipthick
242 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
243 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
244 buflipbot=bordlipbot+lipbufthick
245 bufliptop=bordliptop-lipbufthick
246 if ((lipbufthick*2.0d0).gt.lipthick)
247 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
249 c write(iout,*) "bordliptop=",bordliptop
250 c write(iout,*) "bordlipbot=",bordlipbot
251 c write(iout,*) "bufliptop=",bufliptop
252 c write(iout,*) "buflipbot=",buflipbot
255 if (me.eq.king .or. .not.out1file )
256 & write (iout,*) "DISTCHAINMAX",distchainmax
258 if(me.eq.king.or..not.out1file)
259 & write (iout,'(2a)') diagmeth(kdiag),
260 & ' routine used to diagonalize matrices.'
263 c--------------------------------------------------------------------------
264 subroutine read_REMDpar
268 implicit real*8 (a-h,o-z)
270 include 'COMMON.IOUNITS'
271 include 'COMMON.TIME1'
274 include 'COMMON.LANGEVIN'
276 include 'COMMON.LANGEVIN.lang0'
278 include 'COMMON.INTERACT'
279 include 'COMMON.NAMES'
281 include 'COMMON.REMD'
282 include 'COMMON.CONTROL'
283 include 'COMMON.SETUP'
285 character*320 controlcard
286 character*3200 controlcard1
287 integer iremd_m_total
289 if(me.eq.king.or..not.out1file)
290 & write (iout,*) "REMD setup"
292 call card_concat(controlcard)
293 call readi(controlcard,"NREP",nrep,3)
294 call readi(controlcard,"NSTEX",nstex,1000)
295 call reada(controlcard,"RETMIN",retmin,10.0d0)
296 call reada(controlcard,"RETMAX",retmax,1000.0d0)
297 mremdsync=(index(controlcard,'SYNC').gt.0)
298 call readi(controlcard,"NSYN",i_sync_step,100)
299 restart1file=(index(controlcard,'REST1FILE').gt.0)
300 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
301 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
302 if(max_cache_traj_use.gt.max_cache_traj)
303 & max_cache_traj_use=max_cache_traj
304 if(me.eq.king.or..not.out1file) then
305 cd if (traj1file) then
306 crc caching is in testing - NTWX is not ignored
307 cd write (iout,*) "NTWX value is ignored"
308 cd write (iout,*) " trajectory is stored to one file by master"
309 cd write (iout,*) " before exchange at NSTEX intervals"
311 write (iout,*) "NREP= ",nrep
312 write (iout,*) "NSTEX= ",nstex
313 write (iout,*) "SYNC= ",mremdsync
314 write (iout,*) "NSYN= ",i_sync_step
315 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
318 if (index(controlcard,'TLIST').gt.0) then
320 call card_concat(controlcard1)
321 read(controlcard1,*) (remd_t(i),i=1,nrep)
322 if(me.eq.king.or..not.out1file)
323 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
326 if (index(controlcard,'MLIST').gt.0) then
328 call card_concat(controlcard1)
329 read(controlcard1,*) (remd_m(i),i=1,nrep)
330 if(me.eq.king.or..not.out1file) then
331 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
334 iremd_m_total=iremd_m_total+remd_m(i)
336 write (iout,*) 'Total number of replicas ',iremd_m_total
339 if(me.eq.king.or..not.out1file)
340 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
343 c--------------------------------------------------------------------------
344 subroutine read_MDpar
348 implicit real*8 (a-h,o-z)
350 include 'COMMON.IOUNITS'
351 include 'COMMON.TIME1'
354 include 'COMMON.LANGEVIN'
356 include 'COMMON.LANGEVIN.lang0'
358 include 'COMMON.INTERACT'
359 include 'COMMON.NAMES'
361 include 'COMMON.SETUP'
362 include 'COMMON.CONTROL'
363 include 'COMMON.SPLITELE'
365 character*320 controlcard
367 call card_concat(controlcard)
368 call readi(controlcard,"NSTEP",n_timestep,1000000)
369 call readi(controlcard,"NTWE",ntwe,100)
370 call readi(controlcard,"NTWX",ntwx,1000)
371 call reada(controlcard,"DT",d_time,1.0d-1)
372 call reada(controlcard,"DVMAX",dvmax,2.0d1)
373 call reada(controlcard,"DAMAX",damax,1.0d1)
374 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
375 call readi(controlcard,"LANG",lang,0)
376 RESPA = index(controlcard,"RESPA") .gt. 0
377 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
378 ntime_split0=ntime_split
379 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
380 ntime_split0=ntime_split
381 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
382 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
383 rest = index(controlcard,"REST").gt.0
384 tbf = index(controlcard,"TBF").gt.0
385 usampl = index(controlcard,"USAMPL").gt.0
387 mdpdb = index(controlcard,"MDPDB").gt.0
388 call reada(controlcard,"T_BATH",t_bath,300.0d0)
389 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
390 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
391 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
392 if (count_reset_moment.eq.0) count_reset_moment=1000000000
393 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
394 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
395 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
396 if (count_reset_vel.eq.0) count_reset_vel=1000000000
397 large = index(controlcard,"LARGE").gt.0
398 print_compon = index(controlcard,"PRINT_COMPON").gt.0
399 rattle = index(controlcard,"RATTLE").gt.0
400 preminim = index(controlcard,"PREMINIM").gt.0
402 dccart=(index(controlcard,'CART').gt.0)
405 c if performing umbrella sampling, fragments constrained are read from the fragment file
411 if(me.eq.king.or..not.out1file) then
413 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
415 write (iout,'(a)') "The units are:"
416 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
417 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
418 & " acceleration: angstrom/(48.9 fs)**2"
419 write (iout,'(a)') "energy: kcal/mol, temperature: K"
421 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
422 write (iout,'(a60,f10.5,a)')
423 & "Initial time step of numerical integration:",d_time,
425 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
427 write (iout,'(2a,i4,a)')
428 & "A-MTS algorithm used; initial time step for fast-varying",
429 & " short-range forces split into",ntime_split," steps."
430 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
431 & r_cut," lambda",rlamb
433 write (iout,'(2a,f10.5)')
434 & "Maximum acceleration threshold to reduce the time step",
435 & "/increase split number:",damax
436 write (iout,'(2a,f10.5)')
437 & "Maximum predicted energy drift to reduce the timestep",
438 & "/increase split number:",edriftmax
439 write (iout,'(a60,f10.5)')
440 & "Maximum velocity threshold to reduce velocities:",dvmax
441 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
442 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
443 if (rattle) write (iout,'(a60)')
444 & "Rattle algorithm used to constrain the virtual bonds"
445 if (preminim .or. iranconf.gt.0) then
447 & "Initial structure will be energy-minimized"
452 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
453 call reada(controlcard,"RWAT",rwat,1.4d0)
454 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
455 surfarea=index(controlcard,"SURFAREA").gt.0
456 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
457 if(me.eq.king.or..not.out1file)then
458 write (iout,'(/a,$)') "Langevin dynamics calculation"
461 & " with direct integration of Langevin equations"
462 else if (lang.eq.2) then
463 write (iout,'(a/)') " with TINKER stochasic MD integrator"
464 else if (lang.eq.3) then
465 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
466 else if (lang.eq.4) then
467 write (iout,'(a/)') " in overdamped mode"
469 write (iout,'(//a,i5)')
470 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
473 write (iout,'(a60,f10.5)') "Temperature:",t_bath
474 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
475 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
476 write (iout,'(a60,f10.5)')
477 & "Scaling factor of the friction forces:",scal_fric
478 if (surfarea) write (iout,'(2a,i10,a)')
479 & "Friction coefficients will be scaled by solvent-accessible",
480 & " surface area every",reset_fricmat," steps."
482 c Calculate friction coefficients and bounds of stochastic forces
483 eta=6*pi*cPoise*etawat
484 if(me.eq.king.or..not.out1file)
485 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
487 gamp=scal_fric*(pstok+rwat)*eta
488 stdfp=dsqrt(2*Rb*t_bath/d_time)
490 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
491 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
493 if(me.eq.king.or..not.out1file)then
494 write (iout,'(/2a/)')
495 & "Radii of site types and friction coefficients and std's of",
496 & " stochastic forces of fully exposed sites"
497 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
499 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
500 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
504 if(me.eq.king.or..not.out1file)then
505 write (iout,'(a)') "Berendsen bath calculation"
506 write (iout,'(a60,f10.5)') "Temperature:",t_bath
507 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
509 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
510 & count_reset_moment," steps"
512 & write (iout,'(a,i10,a)')
513 & "Velocities will be reset at random every",count_reset_vel,
517 if(me.eq.king.or..not.out1file)
518 & write (iout,'(a31)') "Microcanonical mode calculation"
520 if(me.eq.king.or..not.out1file)then
521 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
523 write(iout,*) "MD running with constraints."
524 write(iout,*) "Equilibration time ", eq_time, " mtus."
525 write(iout,*) "Constraining ", nfrag," fragments."
526 write(iout,*) "Length of each fragment, weight and q0:"
528 write (iout,*) "Set of restraints #",iset
530 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
531 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
533 write(iout,*) "constraints between ", npair, "fragments."
534 write(iout,*) "constraint pairs, weights and q0:"
536 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
537 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
539 write(iout,*) "angle constraints within ", nfrag_back,
540 & "backbone fragments."
541 write(iout,*) "fragment, weights:"
543 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
544 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
545 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
548 iset=mod(kolor,nset)+1
551 if(me.eq.king.or..not.out1file)
552 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
555 c------------------------------------------------------------------------------
558 C Read molecular data.
560 implicit real*8 (a-h,o-z)
566 include 'COMMON.IOUNITS'
569 include 'COMMON.INTERACT'
570 include 'COMMON.LOCAL'
571 include 'COMMON.NAMES'
572 include 'COMMON.CHAIN'
573 include 'COMMON.FFIELD'
574 include 'COMMON.SBRIDGE'
575 include 'COMMON.HEADER'
576 include 'COMMON.CONTROL'
577 include 'COMMON.DBASE'
578 include 'COMMON.THREAD'
579 include 'COMMON.CONTACTS'
580 include 'COMMON.TORCNSTR'
581 include 'COMMON.TIME1'
582 include 'COMMON.BOUNDS'
584 include 'COMMON.SETUP'
585 character*4 sequence(maxres)
587 double precision x(maxvar)
588 character*256 pdbfile
589 character*320 weightcard
590 character*80 weightcard_t,ucase
591 dimension itype_pdb(maxres)
592 common /pizda/ itype_pdb
593 logical seq_comp,fail
594 double precision energia(0:n_ene)
600 C Read weights of the subsequent energy terms.
601 call card_concat(weightcard)
602 call reada(weightcard,'WLONG',wlong,1.0D0)
603 call reada(weightcard,'WSC',wsc,wlong)
604 call reada(weightcard,'WSCP',wscp,wlong)
605 call reada(weightcard,'WELEC',welec,1.0D0)
606 call reada(weightcard,'WVDWPP',wvdwpp,welec)
607 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
608 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
609 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
610 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
611 call reada(weightcard,'WTURN3',wturn3,1.0D0)
612 call reada(weightcard,'WTURN4',wturn4,1.0D0)
613 call reada(weightcard,'WTURN6',wturn6,1.0D0)
614 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
615 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
616 call reada(weightcard,'WBOND',wbond,1.0D0)
617 call reada(weightcard,'WTOR',wtor,1.0D0)
618 call reada(weightcard,'WTORD',wtor_d,1.0D0)
619 call reada(weightcard,'WANG',wang,1.0D0)
620 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
621 call reada(weightcard,'SCAL14',scal14,0.4D0)
622 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
623 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
624 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
625 call reada(weightcard,'TEMP0',temp0,300.0d0)
626 call reada(weightcard,'WLT',wliptran,0.0D0)
627 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
628 if (index(weightcard,'SOFT').gt.0) ipot=6
629 C 12/1/95 Added weight for the multi-body term WCORR
630 call reada(weightcard,'WCORRH',wcorr,1.0D0)
631 if (wcorr4.gt.0.0d0) wcorr=wcorr4
652 if(me.eq.king.or..not.out1file)
653 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
654 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
656 10 format (/'Energy-term weights (unscaled):'//
657 & 'WSCC= ',f10.6,' (SC-SC)'/
658 & 'WSCP= ',f10.6,' (SC-p)'/
659 & 'WELEC= ',f10.6,' (p-p electr)'/
660 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
661 & 'WBOND= ',f10.6,' (stretching)'/
662 & 'WANG= ',f10.6,' (bending)'/
663 & 'WSCLOC= ',f10.6,' (SC local)'/
664 & 'WTOR= ',f10.6,' (torsional)'/
665 & 'WTORD= ',f10.6,' (double torsional)'/
666 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
667 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
668 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
669 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
670 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
671 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
672 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
673 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
674 & 'WTURN6= ',f10.6,' (turns, 6th order)')
675 if(me.eq.king.or..not.out1file)then
676 if (wcorr4.gt.0.0d0) then
677 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
678 & 'between contact pairs of peptide groups'
679 write (iout,'(2(a,f5.3/))')
680 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
681 & 'Range of quenching the correlation terms:',2*delt_corr
682 else if (wcorr.gt.0.0d0) then
683 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
684 & 'between contact pairs of peptide groups'
686 write (iout,'(a,f8.3)')
687 & 'Scaling factor of 1,4 SC-p interactions:',scal14
688 write (iout,'(a,f8.3)')
689 & 'General scaling factor of SC-p interactions:',scalscp
691 r0_corr=cutoff_corr-delt_corr
693 aad(i,1)=scalscp*aad(i,1)
694 aad(i,2)=scalscp*aad(i,2)
695 bad(i,1)=scalscp*bad(i,1)
696 bad(i,2)=scalscp*bad(i,2)
698 call rescale_weights(t_bath)
699 if(me.eq.king.or..not.out1file)
700 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
701 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
703 22 format (/'Energy-term weights (scaled):'//
704 & 'WSCC= ',f10.6,' (SC-SC)'/
705 & 'WSCP= ',f10.6,' (SC-p)'/
706 & 'WELEC= ',f10.6,' (p-p electr)'/
707 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
708 & 'WBOND= ',f10.6,' (stretching)'/
709 & 'WANG= ',f10.6,' (bending)'/
710 & 'WSCLOC= ',f10.6,' (SC local)'/
711 & 'WTOR= ',f10.6,' (torsional)'/
712 & 'WTORD= ',f10.6,' (double torsional)'/
713 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
714 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
715 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
716 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
717 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
718 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
719 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
720 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
721 & 'WTURN6= ',f10.6,' (turns, 6th order)')
722 if(me.eq.king.or..not.out1file)
723 & write (iout,*) "Reference temperature for weights calculation:",
725 call reada(weightcard,"D0CM",d0cm,3.78d0)
726 call reada(weightcard,"AKCM",akcm,15.1d0)
727 call reada(weightcard,"AKTH",akth,11.0d0)
728 call reada(weightcard,"AKCT",akct,12.0d0)
729 call reada(weightcard,"V1SS",v1ss,-1.08d0)
730 call reada(weightcard,"V2SS",v2ss,7.61d0)
731 call reada(weightcard,"V3SS",v3ss,13.7d0)
732 call reada(weightcard,"EBR",ebr,-5.50D0)
733 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
735 dyn_ss_mask(i)=.false.
739 dyn_ssbond_ij(i,j)=1.0d300
742 call reada(weightcard,"HT",Ht,0.0D0)
744 ss_depth=ebr/wsc-0.25*eps(1,1)
745 Ht=Ht/wsc-0.25*eps(1,1)
746 akcm=akcm*wstrain/wsc
747 akth=akth*wstrain/wsc
748 akct=akct*wstrain/wsc
749 v1ss=v1ss*wstrain/wsc
750 v2ss=v2ss*wstrain/wsc
751 v3ss=v3ss*wstrain/wsc
753 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
756 if(me.eq.king.or..not.out1file) then
757 write (iout,*) "Parameters of the SS-bond potential:"
758 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
760 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
761 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
762 write (iout,*)" HT",Ht
763 print *,'indpdb=',indpdb,' pdbref=',pdbref
765 if (indpdb.gt.0 .or. pdbref) then
766 read(inp,'(a)') pdbfile
767 if(me.eq.king.or..not.out1file)
768 & write (iout,'(2a)') 'PDB data will be read from file ',
769 & pdbfile(:ilen(pdbfile))
770 open(ipdbin,file=pdbfile,status='old',err=33)
772 33 write (iout,'(a)') 'Error opening PDB file.'
775 c print *,'Begin reading pdb data'
784 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
785 & (crefjlee(j,i+nres),j=1,3)
788 c print *,'Finished reading pdb data'
789 if(me.eq.king.or..not.out1file)
790 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
791 & ' nstart_sup=',nstart_sup
793 itype_pdb(i)=itype(i)
797 nct=nstart_sup+nsup-1
798 call contact(.false.,ncont_ref,icont_ref,co)
801 if(me.eq.king.or..not.out1file)
802 & write(iout,*)'Adding sidechains'
806 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
809 do while (fail.and.nsi.le.maxsi)
810 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
814 c Calculalte only the coordinates of the current sidechain; no need to rebuild
816 call locate_side_chain(i)
817 if(fail) write(iout,*)'Adding sidechain failed for res ',
818 & i,' after ',nsi,' trials'
821 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
825 if (indpdb.eq.0) then
826 C Read sequence if not taken from the pdb file.
828 c print *,'nres=',nres
829 if (iscode.gt.0) then
830 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
832 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
834 C Convert sequence to numeric code
836 itype(i)=rescode(i,sequence(i),iscode)
838 C Assign initial virtual bond lengths
844 vbld(i+nres)=dsc(iabs(itype(i)))
845 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
846 c write (iout,*) "i",i," itype",itype(i),
847 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
851 c print '(20i4)',(itype(i),i=1,nres)
854 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
856 if (itype(i).eq.ntyp1) then
860 else if (iabs(itype(i+1)).ne.20) then
862 else if (iabs(itype(i)).ne.20) then
869 if(me.eq.king.or..not.out1file)then
870 write (iout,*) "ITEL"
872 write (iout,*) i,itype(i),itel(i)
874 print *,'Call Read_Bridge.'
877 C 8/13/98 Set limits to generating the dihedral angles
882 read (inp,*) ndih_constr
883 write (iout,*) "ndish_constr",ndih_constr
884 if (ndih_constr.gt.0) then
886 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
887 if(me.eq.king.or..not.out1file)then
889 & 'There are',ndih_constr,' constraints on phi angles.'
891 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
895 phi0(i)=deg2rad*phi0(i)
896 drange(i)=deg2rad*drange(i)
898 if(me.eq.king.or..not.out1file)
899 & write (iout,*) 'FTORS',ftors
902 phibound(1,ii) = phi0(i)-drange(i)
903 phibound(2,ii) = phi0(i)+drange(i)
910 write (iout,'(a)') 'Boundaries in phi angle sampling:'
912 write (iout,'(a3,i5,2f10.1)')
913 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
919 cd print *,'NNT=',NNT,' NCT=',NCT
920 if (itype(1).eq.ntyp1) nnt=2
921 if (itype(nres).eq.ntyp1) nct=nct-1
923 if(me.eq.king.or..not.out1file)
924 & write (iout,'(a,i3)') 'nsup=',nsup
926 if (nsup.le.(nct-nnt+1)) then
927 do i=0,nct-nnt+1-nsup
928 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
934 & 'Error - sequences to be superposed do not match.'
937 do i=0,nsup-(nct-nnt+1)
938 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
940 nstart_sup=nstart_sup+i
946 & 'Error - sequences to be superposed do not match.'
949 if (nsup.eq.0) nsup=nct-nnt
950 if (nstart_sup.eq.0) nstart_sup=nnt
951 if (nstart_seq.eq.0) nstart_seq=nnt
952 if(me.eq.king.or..not.out1file)
953 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
954 & ' nstart_seq=',nstart_seq
956 c--- Zscore rms -------
957 if (nz_start.eq.0) nz_start=nnt
958 if (nz_end.eq.0 .and. nsup.gt.0) then
960 else if (nz_end.eq.0) then
963 if(me.eq.king.or..not.out1file)then
964 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
965 write (iout,*) 'IZ_SC=',iz_sc
967 c----------------------
970 if (.not.pdbref) then
971 call read_angles(inp,*38)
973 38 write (iout,'(a)') 'Error reading reference structure.'
975 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
976 stop 'Error reading reference structure'
978 39 call chainbuild_extconf
980 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
990 call contact(.true.,ncont_ref,icont_ref,co)
992 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
994 if (constr_dist.gt.0) call read_dist_constr
995 c write (iout,*) "After read_dist_constr nhpb",nhpb
996 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
997 if(me.eq.king.or..not.out1file)
998 & write (iout,*) 'Contact order:',co
1000 if(me.eq.king.or..not.out1file)
1001 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1004 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1006 if(me.eq.king.or..not.out1file)
1007 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1008 & icont_ref(1,i),' ',
1009 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1013 write (iout,*) "calling read_saxs_consrtr",nsaxs
1014 if (nsaxs.gt.0) call read_saxs_constr
1017 if (constr_homology.gt.0) then
1018 call read_constr_homology
1019 if (indpdb.gt.0 .or. pdbref) then
1022 c(j,i)=crefjlee(j,i)
1023 cref(j,i,1)=crefjlee(j,i)
1028 write (iout,*) "Array C"
1030 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1031 & (c(j,i+nres),j=1,3)
1033 write (iout,*) "Array Cref"
1035 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
1036 & (cref(j,i+nres,1),j=1,3)
1039 call int_from_cart1(.false.)
1040 call sc_loc_geom(.false.)
1042 thetaref(i)=theta(i)
1047 dc(j,i)=c(j,i+1)-c(j,i)
1048 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1053 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1054 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1062 if (nhpb.gt.0) call hpb_partition
1063 c write (iout,*) "After read_dist_constr nhpb",nhpb
1065 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1066 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1067 & modecalc.ne.10) then
1068 C If input structure hasn't been supplied from the PDB file read or generate
1070 if (iranconf.eq.0 .and. .not. extconf) then
1071 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1072 & write (iout,'(a)') 'Initial geometry will be read in.'
1074 read(inp,'(8f10.5)',end=36,err=36)
1075 & ((c(l,k),l=1,3),k=1,nres),
1076 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1077 write (iout,*) "Exit READ_CART"
1078 write (iout,'(8f10.5)')
1079 & ((c(l,k),l=1,3),k=1,nres),
1080 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1081 call int_from_cart1(.true.)
1082 write (iout,*) "Finish INT_TO_CART"
1085 dc(j,i)=c(j,i+1)-c(j,i)
1086 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1090 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1092 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1093 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1099 write (iout,*) "Calling read_ang"
1100 call read_angles(inp,*36)
1101 write (iout,*) "Calling chainbuild"
1102 call chainbuild_extconf
1105 36 write (iout,'(a)') 'Error reading angle file.'
1107 call mpi_finalize( MPI_COMM_WORLD,IERR )
1109 stop 'Error reading angle file.'
1111 else if (extconf) then
1112 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1113 & write (iout,'(a)') 'Extended chain initial geometry.'
1115 theta(i)=90d0*deg2rad
1118 phi(i)=180d0*deg2rad
1121 alph(i)=110d0*deg2rad
1124 omeg(i)=-120d0*deg2rad
1125 if (itype(i).le.0) omeg(i)=-omeg(i)
1127 c from old chainbuild
1129 C Define the origin and orientation of the coordinate system and locate the
1130 C first three CA's and SC(2).
1134 * Build the alpha-carbon chain.
1137 call locate_next_res(i)
1140 C First and last SC must coincide with the corresponding CA.
1144 dc_norm(j,nres+1)=0.0D0
1145 dc(j,nres+nres)=0.0D0
1146 dc_norm(j,nres+nres)=0.0D0
1148 c(j,nres+nres)=c(j,nres)
1151 C Define the origin and orientation of the coordinate system and locate the
1152 C first three CA's and SC(2).
1156 * Build the alpha-carbon chain.
1159 call locate_next_res(i)
1162 C First and last SC must coincide with the corresponding CA.
1166 dc_norm(j,nres+1)=0.0D0
1167 dc(j,nres+nres)=0.0D0
1168 dc_norm(j,nres+nres)=0.0D0
1170 c(j,nres+nres)=c(j,nres)
1175 if(me.eq.king.or..not.out1file)
1176 & write (iout,'(a)') 'Random-generated initial geometry.'
1180 if (me.eq.king .or. fg_rank.eq.0 .and. (
1181 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1185 call gen_rand_conf(itmp,*30)
1187 30 write (iout,*) 'Failed to generate random conformation',
1188 & ', itrial=',itrial
1189 write (*,*) 'Processor:',me,
1190 & ' Failed to generate random conformation',
1200 write (iout,'(a,i3,a)') 'Processor:',me,
1201 & ' error in generating random conformation.'
1202 write (*,'(a,i3,a)') 'Processor:',me,
1203 & ' error in generating random conformation.'
1206 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1211 & ' error in generating random conformation.'
1216 elseif (modecalc.eq.4) then
1217 read (inp,'(a)') intinname
1218 open (intin,file=intinname,status='old',err=333)
1219 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1220 & write (iout,'(a)') 'intinname',intinname
1221 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1223 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1225 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1227 stop 'Error opening angle file.'
1231 C Generate distance constraints, if the PDB structure is to be regularized.
1232 if (nthread.gt.0) then
1233 call read_threadbase
1236 if (me.eq.king .or. .not. out1file)
1238 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1239 write (iout,'(/a,i3,a)')
1240 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1241 write (iout,'(20i4)') (iss(i),i=1,ns)
1243 write(iout,*)"Running with dynamic disulfide-bond formation"
1245 write (iout,'(/a/)') 'Pre-formed links are:'
1251 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1252 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1258 if (ns.gt.0.and.dyn_ss) then
1262 forcon(i-nss)=forcon(i)
1269 dyn_ss_mask(iss(i))=.true.
1272 if (i2ndstr.gt.0) call secstrp2dihc
1273 c call geom_to_var(nvar,x)
1274 c call etotal(energia(0))
1275 c call enerprint(energia(0))
1276 c call briefout(0,etot)
1278 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1279 cd write (iout,'(a)') 'Variable list:'
1280 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1282 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1283 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1284 & 'Processor',myrank,': end reading molecular data.'
1288 c--------------------------------------------------------------------------
1289 logical function seq_comp(itypea,itypeb,length)
1291 integer length,itypea(length),itypeb(length)
1294 if (itypea(i).ne.itypeb(i)) then
1302 c-----------------------------------------------------------------------------
1303 subroutine read_bridge
1304 C Read information about disulfide bridges.
1305 implicit real*8 (a-h,o-z)
1306 include 'DIMENSIONS'
1310 include 'COMMON.IOUNITS'
1311 include 'COMMON.GEO'
1312 include 'COMMON.VAR'
1313 include 'COMMON.INTERACT'
1314 include 'COMMON.LOCAL'
1315 include 'COMMON.NAMES'
1316 include 'COMMON.CHAIN'
1317 include 'COMMON.FFIELD'
1318 include 'COMMON.SBRIDGE'
1319 include 'COMMON.HEADER'
1320 include 'COMMON.CONTROL'
1321 include 'COMMON.DBASE'
1322 include 'COMMON.THREAD'
1323 include 'COMMON.TIME1'
1324 include 'COMMON.SETUP'
1325 C Read bridging residues.
1326 read (inp,*) ns,(iss(i),i=1,ns)
1328 if(me.eq.king.or..not.out1file)
1329 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1330 C Check whether the specified bridging residues are cystines.
1332 if (itype(iss(i)).ne.1) then
1333 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1334 & 'Do you REALLY think that the residue ',
1335 & restyp(itype(iss(i))),i,
1336 & ' can form a disulfide bridge?!!!'
1337 write (*,'(2a,i3,a)')
1338 & 'Do you REALLY think that the residue ',
1339 & restyp(itype(iss(i))),i,
1340 & ' can form a disulfide bridge?!!!'
1342 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1347 C Read preformed bridges.
1349 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1351 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1354 C Check if the residues involved in bridges are in the specified list of
1355 C bridging residues.
1358 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1359 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1360 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1361 & ' contains residues present in other pairs.'
1362 write (*,'(a,i3,a)') 'Disulfide pair',i,
1363 & ' contains residues present in other pairs.'
1365 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1371 if (ihpb(i).eq.iss(j)) goto 10
1373 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1376 if (jhpb(i).eq.iss(j)) goto 20
1378 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1384 ihpb(i)=ihpb(i)+nres
1385 jhpb(i)=jhpb(i)+nres
1391 c----------------------------------------------------------------------------
1392 subroutine read_x(kanal,*)
1393 implicit real*8 (a-h,o-z)
1394 include 'DIMENSIONS'
1395 include 'COMMON.GEO'
1396 include 'COMMON.VAR'
1397 include 'COMMON.CHAIN'
1398 include 'COMMON.IOUNITS'
1399 include 'COMMON.CONTROL'
1400 include 'COMMON.LOCAL'
1401 include 'COMMON.INTERACT'
1402 c Read coordinates from input
1404 read(kanal,'(8f10.5)',end=10,err=10)
1405 & ((c(l,k),l=1,3),k=1,nres),
1406 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1409 c(j,2*nres)=c(j,nres)
1411 call int_from_cart1(.false.)
1414 dc(j,i)=c(j,i+1)-c(j,i)
1415 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1419 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1421 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1422 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1430 c----------------------------------------------------------------------------
1431 subroutine read_threadbase
1432 implicit real*8 (a-h,o-z)
1433 include 'DIMENSIONS'
1434 include 'COMMON.IOUNITS'
1435 include 'COMMON.GEO'
1436 include 'COMMON.VAR'
1437 include 'COMMON.INTERACT'
1438 include 'COMMON.LOCAL'
1439 include 'COMMON.NAMES'
1440 include 'COMMON.CHAIN'
1441 include 'COMMON.FFIELD'
1442 include 'COMMON.SBRIDGE'
1443 include 'COMMON.HEADER'
1444 include 'COMMON.CONTROL'
1445 include 'COMMON.DBASE'
1446 include 'COMMON.THREAD'
1447 include 'COMMON.TIME1'
1448 C Read pattern database for threading.
1449 read (icbase,*) nseq
1451 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1452 & nres_base(2,i),nres_base(3,i)
1453 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1455 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1456 c & nres_base(2,i),nres_base(3,i)
1457 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1461 if (weidis.eq.0.0D0) weidis=0.1D0
1470 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1471 write (iout,'(a,i5)') 'nexcl: ',nexcl
1472 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1475 c------------------------------------------------------------------------------
1476 subroutine setup_var
1477 implicit real*8 (a-h,o-z)
1478 include 'DIMENSIONS'
1479 include 'COMMON.IOUNITS'
1480 include 'COMMON.GEO'
1481 include 'COMMON.VAR'
1482 include 'COMMON.INTERACT'
1483 include 'COMMON.LOCAL'
1484 include 'COMMON.NAMES'
1485 include 'COMMON.CHAIN'
1486 include 'COMMON.FFIELD'
1487 include 'COMMON.SBRIDGE'
1488 include 'COMMON.HEADER'
1489 include 'COMMON.CONTROL'
1490 include 'COMMON.DBASE'
1491 include 'COMMON.THREAD'
1492 include 'COMMON.TIME1'
1493 C Set up variable list.
1499 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1501 ialph(i,1)=nvar+nside
1505 if (indphi.gt.0) then
1507 else if (indback.gt.0) then
1512 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1515 c----------------------------------------------------------------------------
1516 subroutine gen_dist_constr
1517 C Generate CA distance constraints.
1518 implicit real*8 (a-h,o-z)
1519 include 'DIMENSIONS'
1520 include 'COMMON.IOUNITS'
1521 include 'COMMON.GEO'
1522 include 'COMMON.VAR'
1523 include 'COMMON.INTERACT'
1524 include 'COMMON.LOCAL'
1525 include 'COMMON.NAMES'
1526 include 'COMMON.CHAIN'
1527 include 'COMMON.FFIELD'
1528 include 'COMMON.SBRIDGE'
1529 include 'COMMON.HEADER'
1530 include 'COMMON.CONTROL'
1531 include 'COMMON.DBASE'
1532 include 'COMMON.THREAD'
1533 include 'COMMON.TIME1'
1534 dimension itype_pdb(maxres)
1535 common /pizda/ itype_pdb
1537 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1538 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1539 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1541 do i=nstart_sup,nstart_sup+nsup-1
1542 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1543 cd & ' seq_pdb', restyp(itype_pdb(i))
1544 do j=i+2,nstart_sup+nsup-1
1546 ihpb(nhpb)=i+nstart_seq-nstart_sup
1547 jhpb(nhpb)=j+nstart_seq-nstart_sup
1549 dhpb(nhpb)=dist(i,j)
1552 cd write (iout,'(a)') 'Distance constraints:'
1557 cd if (ii.gt.nres) then
1562 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1563 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1564 cd & dhpb(i),forcon(i)
1568 c----------------------------------------------------------------------------
1570 implicit real*8 (a-h,o-z)
1571 include 'DIMENSIONS'
1572 include 'COMMON.MAP'
1573 include 'COMMON.IOUNITS'
1574 character*3 angid(4) /'THE','PHI','ALP','OME'/
1575 character*80 mapcard,ucase
1577 read (inp,'(a)') mapcard
1578 mapcard=ucase(mapcard)
1579 if (index(mapcard,'PHI').gt.0) then
1581 else if (index(mapcard,'THE').gt.0) then
1583 else if (index(mapcard,'ALP').gt.0) then
1585 else if (index(mapcard,'OME').gt.0) then
1588 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1589 stop 'Error - illegal variable spec in MAP card.'
1591 call readi (mapcard,'RES1',res1(imap),0)
1592 call readi (mapcard,'RES2',res2(imap),0)
1593 if (res1(imap).eq.0) then
1594 res1(imap)=res2(imap)
1595 else if (res2(imap).eq.0) then
1596 res2(imap)=res1(imap)
1598 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1600 & 'Error - illegal definition of variable group in MAP.'
1601 stop 'Error - illegal definition of variable group in MAP.'
1603 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1604 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1605 call readi(mapcard,'NSTEP',nstep(imap),0)
1606 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1608 & 'Illegal boundary and/or step size specification in MAP.'
1609 stop 'Illegal boundary and/or step size specification in MAP.'
1614 c----------------------------------------------------------------------------
1616 implicit real*8 (a-h,o-z)
1617 include 'DIMENSIONS'
1618 include 'COMMON.IOUNITS'
1619 include 'COMMON.GEO'
1620 include 'COMMON.CSA'
1621 include 'COMMON.BANK'
1622 include 'COMMON.CONTROL'
1624 character*620 mcmcard
1625 call card_concat(mcmcard)
1627 call readi(mcmcard,'NCONF',nconf,50)
1628 call readi(mcmcard,'NADD',nadd,0)
1629 call readi(mcmcard,'JSTART',jstart,1)
1630 call readi(mcmcard,'JEND',jend,1)
1631 call readi(mcmcard,'NSTMAX',nstmax,500000)
1632 call readi(mcmcard,'N0',n0,1)
1633 call readi(mcmcard,'N1',n1,6)
1634 call readi(mcmcard,'N2',n2,4)
1635 call readi(mcmcard,'N3',n3,0)
1636 call readi(mcmcard,'N4',n4,0)
1637 call readi(mcmcard,'N5',n5,0)
1638 call readi(mcmcard,'N6',n6,10)
1639 call readi(mcmcard,'N7',n7,0)
1640 call readi(mcmcard,'N8',n8,0)
1641 call readi(mcmcard,'N9',n9,0)
1642 call readi(mcmcard,'N14',n14,0)
1643 call readi(mcmcard,'N15',n15,0)
1644 call readi(mcmcard,'N16',n16,0)
1645 call readi(mcmcard,'N17',n17,0)
1646 call readi(mcmcard,'N18',n18,0)
1648 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1650 call readi(mcmcard,'NDIFF',ndiff,2)
1651 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1652 call readi(mcmcard,'IS1',is1,1)
1653 call readi(mcmcard,'IS2',is2,8)
1654 call readi(mcmcard,'NRAN0',nran0,4)
1655 call readi(mcmcard,'NRAN1',nran1,2)
1656 call readi(mcmcard,'IRR',irr,1)
1657 call readi(mcmcard,'NSEED',nseed,20)
1658 call readi(mcmcard,'NTOTAL',ntotal,10000)
1659 call reada(mcmcard,'CUT1',cut1,2.0d0)
1660 call reada(mcmcard,'CUT2',cut2,5.0d0)
1661 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1662 call readi(mcmcard,'ICMAX',icmax,3)
1663 call readi(mcmcard,'IRESTART',irestart,0)
1664 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1667 call reada(mcmcard,'DELE',dele,20.0d0)
1668 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1669 call readi(mcmcard,'IREF',iref,0)
1670 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1671 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1672 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1673 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1674 write (iout,*) "NCONF_IN",nconf_in
1677 c----------------------------------------------------------------------------
1678 cfmc subroutine mcmfread
1679 cfmc implicit real*8 (a-h,o-z)
1680 cfmc include 'DIMENSIONS'
1681 cfmc include 'COMMON.MCMF'
1682 cfmc include 'COMMON.IOUNITS'
1683 cfmc include 'COMMON.GEO'
1684 cfmc character*80 ucase
1685 cfmc character*620 mcmcard
1686 cfmc call card_concat(mcmcard)
1688 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1689 cfmc write(iout,*)'MAXRANT=',maxrant
1690 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1691 cfmc write(iout,*)'MAXFAM=',maxfam
1692 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1693 cfmc write(iout,*)'NNET1=',nnet1
1694 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1695 cfmc write(iout,*)'NNET2=',nnet2
1696 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1697 cfmc write(iout,*)'NNET3=',nnet3
1698 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1699 cfmc write(iout,*)'ILASTT=',ilastt
1700 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1701 cfmc write(iout,*)'MAXSTR=',maxstr
1702 cfmc maxstr_f=maxstr/maxfam
1703 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1704 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1705 cfmc write(iout,*)'NMCMF=',nmcmf
1706 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1707 cfmc write(iout,*)'IFOCUS=',ifocus
1708 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1709 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1710 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1711 cfmc write(iout,*)'INTPRT=',intprt
1712 cfmc call readi(mcmcard,'IPRT',iprt,100)
1713 cfmc write(iout,*)'IPRT=',iprt
1714 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1715 cfmc write(iout,*)'IMAXTR=',imaxtr
1716 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1717 cfmc write(iout,*)'MAXEVEN=',maxeven
1718 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1719 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1720 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1721 cfmc write(iout,*)'INIMIN=',inimin
1722 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1723 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1724 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1725 cfmc write(iout,*)'NTHREAD=',nthread
1726 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1727 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1728 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1729 cfmc write(iout,*)'MAXPERT=',maxpert
1730 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1731 cfmc write(iout,*)'IRMSD=',irmsd
1732 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1733 cfmc write(iout,*)'DENEMIN=',denemin
1734 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1735 cfmc write(iout,*)'RCUT1S=',rcut1s
1736 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1737 cfmc write(iout,*)'RCUT1E=',rcut1e
1738 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1739 cfmc write(iout,*)'RCUT2S=',rcut2s
1740 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1741 cfmc write(iout,*)'RCUT2E=',rcut2e
1742 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1743 cfmc write(iout,*)'DPERT1=',d_pert1
1744 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1745 cfmc write(iout,*)'DPERT1A=',d_pert1a
1746 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1747 cfmc write(iout,*)'DPERT2=',d_pert2
1748 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1749 cfmc write(iout,*)'DPERT2A=',d_pert2a
1750 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1751 cfmc write(iout,*)'DPERT2B=',d_pert2b
1752 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1753 cfmc write(iout,*)'DPERT2C=',d_pert2c
1754 cfmc d_pert1=deg2rad*d_pert1
1755 cfmc d_pert1a=deg2rad*d_pert1a
1756 cfmc d_pert2=deg2rad*d_pert2
1757 cfmc d_pert2a=deg2rad*d_pert2a
1758 cfmc d_pert2b=deg2rad*d_pert2b
1759 cfmc d_pert2c=deg2rad*d_pert2c
1760 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1761 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1762 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1763 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1764 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1765 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1766 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1767 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1768 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1769 cfmc write(iout,*)'RCUTINI=',rcutini
1770 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1771 cfmc write(iout,*)'GRAT=',grat
1772 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1773 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1777 c----------------------------------------------------------------------------
1779 implicit real*8 (a-h,o-z)
1780 include 'DIMENSIONS'
1781 include 'COMMON.MCM'
1782 include 'COMMON.MCE'
1783 include 'COMMON.IOUNITS'
1785 character*320 mcmcard
1786 call card_concat(mcmcard)
1787 call readi(mcmcard,'MAXACC',maxacc,100)
1788 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1789 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1790 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1791 call readi(mcmcard,'MAXREPM',maxrepm,200)
1792 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1793 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1794 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1795 call reada(mcmcard,'E_UP',e_up,5.0D0)
1796 call reada(mcmcard,'DELTE',delte,0.1D0)
1797 call readi(mcmcard,'NSWEEP',nsweep,5)
1798 call readi(mcmcard,'NSTEPH',nsteph,0)
1799 call readi(mcmcard,'NSTEPC',nstepc,0)
1800 call reada(mcmcard,'TMIN',tmin,298.0D0)
1801 call reada(mcmcard,'TMAX',tmax,298.0D0)
1802 call readi(mcmcard,'NWINDOW',nwindow,0)
1803 call readi(mcmcard,'PRINT_MC',print_mc,0)
1804 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1805 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1806 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1807 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1808 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1809 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1810 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1811 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1812 if (nwindow.gt.0) then
1813 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1815 winlen(i)=winend(i)-winstart(i)+1
1818 if (tmax.lt.tmin) tmax=tmin
1819 if (tmax.eq.tmin) then
1823 if (nstepc.gt.0 .and. nsteph.gt.0) then
1824 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1825 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1827 C Probabilities of different move types
1828 sumpro_type(0)=0.0D0
1829 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1830 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1831 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1832 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1833 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1834 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1835 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1837 print *,'i',i,' sumprotype',sumpro_type(i)
1838 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1839 print *,'i',i,' sumprotype',sumpro_type(i)
1843 c----------------------------------------------------------------------------
1844 subroutine read_minim
1845 implicit real*8 (a-h,o-z)
1846 include 'DIMENSIONS'
1847 include 'COMMON.MINIM'
1848 include 'COMMON.IOUNITS'
1849 include 'COMMON.CONTROL'
1850 include 'COMMON.SETUP'
1852 character*320 minimcard
1853 call card_concat(minimcard)
1854 call readi(minimcard,'MAXMIN',maxmin,2000)
1855 call readi(minimcard,'MAXFUN',maxfun,5000)
1856 call readi(minimcard,'MINMIN',minmin,maxmin)
1857 call readi(minimcard,'MINFUN',minfun,maxmin)
1858 call reada(minimcard,'TOLF',tolf,1.0D-2)
1859 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1860 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1861 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1862 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1864 if (.not. out1file .or. me.eq.king) then
1866 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1867 & 'Options in energy minimization:'
1868 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1869 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1870 & 'MinMin:',MinMin,' MinFun:',MinFun,
1871 & ' TolF:',TolF,' RTolF:',RTolF
1877 c----------------------------------------------------------------------------
1878 subroutine read_angles(kanal,*)
1879 implicit real*8 (a-h,o-z)
1880 include 'DIMENSIONS'
1881 include 'COMMON.GEO'
1882 include 'COMMON.VAR'
1883 include 'COMMON.CHAIN'
1884 include 'COMMON.IOUNITS'
1885 include 'COMMON.CONTROL'
1886 c Read angles from input
1888 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1889 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1890 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1891 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1894 c 9/7/01 avoid 180 deg valence angle
1895 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1897 theta(i)=deg2rad*theta(i)
1898 phi(i)=deg2rad*phi(i)
1899 alph(i)=deg2rad*alph(i)
1900 omeg(i)=deg2rad*omeg(i)
1905 c----------------------------------------------------------------------------
1906 subroutine reada(rekord,lancuch,wartosc,default)
1908 character*(*) rekord,lancuch
1909 double precision wartosc,default
1912 iread=index(rekord,lancuch)
1913 if (iread.eq.0) then
1917 iread=iread+ilen(lancuch)+1
1918 read (rekord(iread:),*,err=10,end=10) wartosc
1923 c----------------------------------------------------------------------------
1924 subroutine readi(rekord,lancuch,wartosc,default)
1926 character*(*) rekord,lancuch
1927 integer wartosc,default
1930 iread=index(rekord,lancuch)
1931 if (iread.eq.0) then
1935 iread=iread+ilen(lancuch)+1
1936 read (rekord(iread:),*,err=10,end=10) wartosc
1941 c----------------------------------------------------------------------------
1942 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1945 integer tablica(dim),default
1946 character*(*) rekord,lancuch
1953 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1954 if (iread.eq.0) return
1955 iread=iread+ilen(lancuch)+1
1956 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1959 c----------------------------------------------------------------------------
1960 subroutine multreada(rekord,lancuch,tablica,dim,default)
1963 double precision tablica(dim),default
1964 character*(*) rekord,lancuch
1971 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1972 if (iread.eq.0) return
1973 iread=iread+ilen(lancuch)+1
1974 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1977 c----------------------------------------------------------------------------
1978 subroutine openunits
1979 implicit real*8 (a-h,o-z)
1980 include 'DIMENSIONS'
1983 character*16 form,nodename
1986 include 'COMMON.SETUP'
1987 include 'COMMON.IOUNITS'
1989 include 'COMMON.CONTROL'
1990 integer lenpre,lenpot,ilen,lentmp
1992 character*3 out1file_text,ucase
1995 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1996 call getenv_loc("PREFIX",prefix)
1998 call getenv_loc("POT",pot)
1999 call getenv_loc("DIRTMP",tmpdir)
2000 call getenv_loc("CURDIR",curdir)
2001 call getenv_loc("OUT1FILE",out1file_text)
2002 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2003 out1file_text=ucase(out1file_text)
2004 if (out1file_text(1:1).eq."Y") then
2007 out1file=fg_rank.gt.0
2012 if (lentmp.gt.0) then
2013 write (*,'(80(1h!))')
2014 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2015 write (*,'(80(1h!))')
2016 write (*,*)"All output files will be on node /tmp directory."
2018 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2019 if (me.eq.king) then
2020 write (*,*) "The master node is ",nodename
2021 else if (fg_rank.eq.0) then
2022 write (*,*) "I am the CG slave node ",nodename
2024 write (*,*) "I am the FG slave node ",nodename
2027 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2028 lenpre = lentmp+lenpre+1
2030 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2031 C Get the names and open the input files
2032 #if defined(WINIFL) || defined(WINPGI)
2033 open(1,file=pref_orig(:ilen(pref_orig))//
2034 & '.inp',status='old',readonly,shared)
2035 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2036 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2037 C Get parameter filenames and open the parameter files.
2038 call getenv_loc('BONDPAR',bondname)
2039 open (ibond,file=bondname,status='old',readonly,shared)
2040 call getenv_loc('THETPAR',thetname)
2041 open (ithep,file=thetname,status='old',readonly,shared)
2043 call getenv_loc('THETPARPDB',thetname_pdb)
2044 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2046 call getenv_loc('ROTPAR',rotname)
2047 open (irotam,file=rotname,status='old',readonly,shared)
2049 call getenv_loc('ROTPARPDB',rotname_pdb)
2050 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2052 call getenv_loc('TORPAR',torname)
2053 open (itorp,file=torname,status='old',readonly,shared)
2054 call getenv_loc('TORDPAR',tordname)
2055 open (itordp,file=tordname,status='old',readonly,shared)
2056 call getenv_loc('FOURIER',fouriername)
2057 open (ifourier,file=fouriername,status='old',readonly,shared)
2058 call getenv_loc('ELEPAR',elename)
2059 open (ielep,file=elename,status='old',readonly,shared)
2060 call getenv_loc('SIDEPAR',sidename)
2061 open (isidep,file=sidename,status='old',readonly,shared)
2062 #elif (defined CRAY) || (defined AIX)
2063 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2065 c print *,"Processor",myrank," opened file 1"
2066 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2067 c print *,"Processor",myrank," opened file 9"
2068 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2069 C Get parameter filenames and open the parameter files.
2070 call getenv_loc('BONDPAR',bondname)
2071 open (ibond,file=bondname,status='old',action='read')
2072 c print *,"Processor",myrank," opened file IBOND"
2073 call getenv_loc('THETPAR',thetname)
2074 open (ithep,file=thetname,status='old',action='read')
2075 c print *,"Processor",myrank," opened file ITHEP"
2077 call getenv_loc('THETPARPDB',thetname_pdb)
2078 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2080 call getenv_loc('ROTPAR',rotname)
2081 open (irotam,file=rotname,status='old',action='read')
2082 c print *,"Processor",myrank," opened file IROTAM"
2084 call getenv_loc('ROTPARPDB',rotname_pdb)
2085 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2087 call getenv_loc('TORPAR',torname)
2088 open (itorp,file=torname,status='old',action='read')
2089 c print *,"Processor",myrank," opened file ITORP"
2090 call getenv_loc('TORDPAR',tordname)
2091 open (itordp,file=tordname,status='old',action='read')
2092 c print *,"Processor",myrank," opened file ITORDP"
2093 call getenv_loc('SCCORPAR',sccorname)
2094 open (isccor,file=sccorname,status='old',action='read')
2095 c print *,"Processor",myrank," opened file ISCCOR"
2096 call getenv_loc('FOURIER',fouriername)
2097 open (ifourier,file=fouriername,status='old',action='read')
2098 c print *,"Processor",myrank," opened file IFOURIER"
2099 call getenv_loc('ELEPAR',elename)
2100 open (ielep,file=elename,status='old',action='read')
2101 c print *,"Processor",myrank," opened file IELEP"
2102 call getenv_loc('SIDEPAR',sidename)
2103 open (isidep,file=sidename,status='old',action='read')
2104 c print *,"Processor",myrank," opened file ISIDEP"
2105 c print *,"Processor",myrank," opened parameter files"
2107 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2108 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2109 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2110 C Get parameter filenames and open the parameter files.
2111 call getenv_loc('BONDPAR',bondname)
2112 open (ibond,file=bondname,status='old')
2113 call getenv_loc('THETPAR',thetname)
2114 open (ithep,file=thetname,status='old')
2116 call getenv_loc('THETPARPDB',thetname_pdb)
2117 open (ithep_pdb,file=thetname_pdb,status='old')
2119 call getenv_loc('ROTPAR',rotname)
2120 open (irotam,file=rotname,status='old')
2122 call getenv_loc('ROTPARPDB',rotname_pdb)
2123 open (irotam_pdb,file=rotname_pdb,status='old')
2125 call getenv_loc('TORPAR',torname)
2126 open (itorp,file=torname,status='old')
2127 call getenv_loc('TORDPAR',tordname)
2128 open (itordp,file=tordname,status='old')
2129 call getenv_loc('SCCORPAR',sccorname)
2130 open (isccor,file=sccorname,status='old')
2131 call getenv_loc('FOURIER',fouriername)
2132 open (ifourier,file=fouriername,status='old')
2133 call getenv_loc('ELEPAR',elename)
2134 open (ielep,file=elename,status='old')
2135 call getenv_loc('SIDEPAR',sidename)
2136 open (isidep,file=sidename,status='old')
2137 call getenv_loc('LIPTRANPAR',liptranname)
2138 open (iliptranpar,file=liptranname,status='old')
2140 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2142 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2143 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2144 C Get parameter filenames and open the parameter files.
2145 call getenv_loc('BONDPAR',bondname)
2146 open (ibond,file=bondname,status='old',readonly)
2147 call getenv_loc('THETPAR',thetname)
2148 open (ithep,file=thetname,status='old',readonly)
2149 call getenv_loc('ROTPAR',rotname)
2150 open (irotam,file=rotname,status='old',readonly)
2151 call getenv_loc('TORPAR',torname)
2152 open (itorp,file=torname,status='old',readonly)
2153 call getenv_loc('TORDPAR',tordname)
2154 open (itordp,file=tordname,status='old',readonly)
2155 call getenv_loc('SCCORPAR',sccorname)
2156 open (isccor,file=sccorname,status='old',readonly)
2158 call getenv_loc('THETPARPDB',thetname_pdb)
2159 c print *,"thetname_pdb ",thetname_pdb
2160 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2161 c print *,ithep_pdb," opened"
2163 call getenv_loc('FOURIER',fouriername)
2164 open (ifourier,file=fouriername,status='old',readonly)
2165 call getenv_loc('ELEPAR',elename)
2166 open (ielep,file=elename,status='old',readonly)
2167 call getenv_loc('SIDEPAR',sidename)
2168 open (isidep,file=sidename,status='old',readonly)
2169 call getenv_loc('LIPTRANPAR',liptranname)
2170 open (iliptranpar,file=liptranname,status='old',readonly)
2172 call getenv_loc('ROTPARPDB',rotname_pdb)
2173 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2178 C 8/9/01 In the newest version SCp interaction constants are read from a file
2179 C Use -DOLDSCP to use hard-coded constants instead.
2181 call getenv_loc('SCPPAR',scpname)
2182 #if defined(WINIFL) || defined(WINPGI)
2183 open (iscpp,file=scpname,status='old',readonly,shared)
2184 #elif (defined CRAY) || (defined AIX)
2185 open (iscpp,file=scpname,status='old',action='read')
2187 open (iscpp,file=scpname,status='old')
2189 open (iscpp,file=scpname,status='old',readonly)
2192 call getenv_loc('PATTERN',patname)
2193 #if defined(WINIFL) || defined(WINPGI)
2194 open (icbase,file=patname,status='old',readonly,shared)
2195 #elif (defined CRAY) || (defined AIX)
2196 open (icbase,file=patname,status='old',action='read')
2198 open (icbase,file=patname,status='old')
2200 open (icbase,file=patname,status='old',readonly)
2203 C Open output file only for CG processes
2204 c print *,"Processor",myrank," fg_rank",fg_rank
2205 if (fg_rank.eq.0) then
2207 if (nodes.eq.1) then
2210 npos = dlog10(dfloat(nodes-1))+1
2212 if (npos.lt.3) npos=3
2213 write (liczba,'(i1)') npos
2214 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2216 write (liczba,form) me
2217 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2218 & liczba(:ilen(liczba))
2219 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2221 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2223 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2224 & liczba(:ilen(liczba))//'.mol2'
2225 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2226 & liczba(:ilen(liczba))//'.stat'
2228 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2229 & //liczba(:ilen(liczba))//'.stat')
2230 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2233 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2234 & liczba(:ilen(liczba))//'.const'
2239 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2240 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2241 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2242 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2243 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2245 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2247 rest2name=prefix(:ilen(prefix))//'.rst'
2249 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2252 #if defined(AIX) || defined(PGI)
2253 if (me.eq.king .or. .not. out1file)
2254 & open(iout,file=outname,status='unknown')
2256 if (fg_rank.gt.0) then
2257 write (liczba,'(i3.3)') myrank/nfgtasks
2258 write (ll,'(bz,i3.3)') fg_rank
2259 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2264 open(igeom,file=intname,status='unknown',position='append')
2265 open(ipdb,file=pdbname,status='unknown')
2266 open(imol2,file=mol2name,status='unknown')
2267 open(istat,file=statname,status='unknown',position='append')
2269 c1out open(iout,file=outname,status='unknown')
2272 if (me.eq.king .or. .not.out1file)
2273 & open(iout,file=outname,status='unknown')
2275 if (fg_rank.gt.0) then
2276 write (liczba,'(i3.3)') myrank/nfgtasks
2277 write (ll,'(bz,i3.3)') fg_rank
2278 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2283 open(igeom,file=intname,status='unknown',access='append')
2284 open(ipdb,file=pdbname,status='unknown')
2285 open(imol2,file=mol2name,status='unknown')
2286 open(istat,file=statname,status='unknown',access='append')
2288 c1out open(iout,file=outname,status='unknown')
2291 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2292 csa_seed=prefix(:lenpre)//'.CSA.seed'
2293 csa_history=prefix(:lenpre)//'.CSA.history'
2294 csa_bank=prefix(:lenpre)//'.CSA.bank'
2295 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2296 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2297 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2298 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2299 csa_int=prefix(:lenpre)//'.int'
2300 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2301 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2302 csa_in=prefix(:lenpre)//'.CSA.in'
2303 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2306 write (iout,'(80(1h-))')
2307 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2308 write (iout,'(80(1h-))')
2309 write (iout,*) "Input file : ",
2310 & pref_orig(:ilen(pref_orig))//'.inp'
2311 write (iout,*) "Output file : ",
2312 & outname(:ilen(outname))
2314 write (iout,*) "Sidechain potential file : ",
2315 & sidename(:ilen(sidename))
2317 write (iout,*) "SCp potential file : ",
2318 & scpname(:ilen(scpname))
2320 write (iout,*) "Electrostatic potential file : ",
2321 & elename(:ilen(elename))
2322 write (iout,*) "Cumulant coefficient file : ",
2323 & fouriername(:ilen(fouriername))
2324 write (iout,*) "Torsional parameter file : ",
2325 & torname(:ilen(torname))
2326 write (iout,*) "Double torsional parameter file : ",
2327 & tordname(:ilen(tordname))
2328 write (iout,*) "SCCOR parameter file : ",
2329 & sccorname(:ilen(sccorname))
2330 write (iout,*) "Bond & inertia constant file : ",
2331 & bondname(:ilen(bondname))
2332 write (iout,*) "Bending parameter file : ",
2333 & thetname(:ilen(thetname))
2334 write (iout,*) "Rotamer parameter file : ",
2335 & rotname(:ilen(rotname))
2336 write (iout,*) "Thetpdb parameter file : ",
2337 & thetname_pdb(:ilen(thetname_pdb))
2338 write (iout,*) "Threading database : ",
2339 & patname(:ilen(patname))
2341 &write (iout,*)" DIRTMP : ",
2343 write (iout,'(80(1h-))')
2347 c----------------------------------------------------------------------------
2348 subroutine card_concat(card)
2349 implicit real*8 (a-h,o-z)
2350 include 'DIMENSIONS'
2351 include 'COMMON.IOUNITS'
2353 character*80 karta,ucase
2355 read (inp,'(a)') karta
2358 do while (karta(80:80).eq.'&')
2359 card=card(:ilen(card)+1)//karta(:79)
2360 read (inp,'(a)') karta
2363 card=card(:ilen(card)+1)//karta
2366 c----------------------------------------------------------------------------------
2368 implicit real*8 (a-h,o-z)
2369 include 'DIMENSIONS'
2370 include 'COMMON.CHAIN'
2371 include 'COMMON.IOUNITS'
2373 include 'COMMON.CONTROL'
2374 open(irest2,file=rest2name,status='unknown')
2375 read(irest2,*) totT,EK,potE,totE,t_bath
2378 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2381 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2383 if(usampl.or.homol_nset.gt.1) then
2384 read (irest2,*) iset
2389 c---------------------------------------------------------------------------------
2390 subroutine read_fragments
2391 implicit real*8 (a-h,o-z)
2392 include 'DIMENSIONS'
2396 include 'COMMON.SETUP'
2397 include 'COMMON.CHAIN'
2398 include 'COMMON.IOUNITS'
2400 include 'COMMON.CONTROL'
2402 read(inp,*) nset,nfrag,npair,nfrag_back
2403 if(me.eq.king.or..not.out1file)
2404 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2405 & " nfrag_back",nfrag_back
2407 read(inp,*) mset(iset1)
2409 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2411 if(me.eq.king.or..not.out1file)
2412 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2413 & ifrag(2,i,iset1), qinfrag(i,iset1)
2416 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2418 if(me.eq.king.or..not.out1file)
2419 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2420 & ipair(2,i,iset1), qinpair(i,iset1)
2423 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2424 & wfrag_back(3,i,iset1),
2425 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2426 if(me.eq.king.or..not.out1file)
2427 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2428 & wfrag_back(2,i,iset1),
2429 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2430 & ifrag_back(2,i,iset1)
2435 C---------------------------------------------------------------------------
2436 subroutine read_afminp
2437 implicit real*8 (a-h,o-z)
2438 include 'DIMENSIONS'
2442 include 'COMMON.SETUP'
2443 include 'COMMON.CONTROL'
2444 include 'COMMON.CHAIN'
2445 include 'COMMON.IOUNITS'
2446 include 'COMMON.SBRIDGE'
2447 character*320 afmcard
2448 c print *, "wchodze"
2449 call card_concat(afmcard)
2450 call readi(afmcard,"BEG",afmbeg,0)
2451 call readi(afmcard,"END",afmend,0)
2452 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2453 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2454 print *,'FORCE=' ,forceAFMconst
2455 CCCC NOW PROPERTIES FOR AFM
2458 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2460 distafminit=dsqrt(distafminit)
2461 print *,'initdist',distafminit
2465 c-------------------------------------------------------------------------------
2466 subroutine read_saxs_constr
2467 implicit real*8 (a-h,o-z)
2468 include 'DIMENSIONS'
2472 include 'COMMON.SETUP'
2473 include 'COMMON.CONTROL'
2474 include 'COMMON.CHAIN'
2475 include 'COMMON.IOUNITS'
2476 include 'COMMON.SBRIDGE'
2477 double precision cm(3)
2479 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2481 if (saxs_mode.eq.0) then
2482 c SAXS distance distribution
2484 read(inp,*) distsaxs(i),Psaxs(i)
2488 Cnorm = Cnorm + Psaxs(i)
2490 write (iout,*) "Cnorm",Cnorm
2492 Psaxs(i)=Psaxs(i)/Cnorm
2494 write (iout,*) "Normalized distance distribution from SAXS"
2496 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2501 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2508 cm(j)=cm(j)+Csaxs(j,i)
2516 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2519 write (iout,*) "SAXS sphere coordinates"
2521 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2526 c-------------------------------------------------------------------------------
2527 subroutine read_dist_constr
2528 implicit real*8 (a-h,o-z)
2529 include 'DIMENSIONS'
2533 include 'COMMON.SETUP'
2534 include 'COMMON.CONTROL'
2535 include 'COMMON.CHAIN'
2536 include 'COMMON.IOUNITS'
2537 include 'COMMON.SBRIDGE'
2538 integer ifrag_(2,100),ipair_(2,100)
2539 double precision wfrag_(100),wpair_(100)
2540 character*500 controlcard
2541 c write (iout,*) "Calling read_dist_constr"
2542 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2544 call card_concat(controlcard)
2545 call readi(controlcard,"NFRAG",nfrag_,0)
2546 call readi(controlcard,"NPAIR",npair_,0)
2547 call readi(controlcard,"NDIST",ndist_,0)
2548 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2549 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2550 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2551 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2552 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2553 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2554 c write (iout,*) "IFRAG"
2556 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2558 c write (iout,*) "IPAIR"
2560 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2562 if (.not.refstr .and. nfrag.gt.0) then
2564 & "ERROR: no reference structure to compute distance restraints"
2566 & "Restraints must be specified explicitly (NDIST=number)"
2569 if (nfrag.lt.2 .and. npair.gt.0) then
2570 write (iout,*) "ERROR: Less than 2 fragments specified",
2571 & " but distance restraints between pairs requested"
2576 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2577 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2578 & ifrag_(2,i)=nstart_sup+nsup-1
2579 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2581 if (wfrag_(i).gt.0.0d0) then
2582 do j=ifrag_(1,i),ifrag_(2,i)-1
2583 do k=j+1,ifrag_(2,i)
2584 c write (iout,*) "j",j," k",k
2586 if (constr_dist.eq.1) then
2591 forcon(nhpb)=wfrag_(i)
2592 else if (constr_dist.eq.2) then
2593 if (ddjk.le.dist_cut) then
2598 forcon(nhpb)=wfrag_(i)
2605 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2608 if (.not.out1file .or. me.eq.king)
2609 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2610 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2612 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2613 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2620 if (wpair_(i).gt.0.0d0) then
2628 do j=ifrag_(1,ii),ifrag_(2,ii)
2629 do k=ifrag_(1,jj),ifrag_(2,jj)
2633 forcon(nhpb)=wpair_(i)
2634 dhpb(nhpb)=dist(j,k)
2636 if (.not.out1file .or. me.eq.king)
2637 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2638 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2640 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2641 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2648 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2649 if (forcon(nhpb+1).gt.0.0d0) then
2651 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2653 if (.not.out1file .or. me.eq.king)
2654 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2655 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2657 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2658 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2665 c-------------------------------------------------------------------------------
2667 subroutine read_constr_homology
2669 include 'DIMENSIONS'
2673 include 'COMMON.SETUP'
2674 include 'COMMON.CONTROL'
2675 include 'COMMON.CHAIN'
2676 include 'COMMON.IOUNITS'
2678 include 'COMMON.GEO'
2679 include 'COMMON.INTERACT'
2680 include 'COMMON.NAMES'
2682 c For new homol impl
2684 include 'COMMON.VAR'
2687 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2689 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2690 c & sigma_odl_temp(maxres,maxres,max_template)
2692 character*24 model_ki_dist, model_ki_angle
2693 character*500 controlcard
2694 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2695 logical lprn /.true./
2700 c FP - Nov. 2014 Temporary specifications for new vars
2702 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2704 double precision, dimension (max_template,maxres) :: rescore
2705 double precision, dimension (max_template,maxres) :: rescore2
2706 double precision, dimension (max_template,maxres) :: rescore3
2707 character*24 pdbfile,tpl_k_rescore
2708 c -----------------------------------------------------------------
2709 c Reading multiple PDB ref structures and calculation of retraints
2710 c not using pre-computed ones stored in files model_ki_{dist,angle}
2712 c -----------------------------------------------------------------
2715 c Alternative: reading from input
2716 call card_concat(controlcard)
2717 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2718 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2719 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2720 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2721 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2722 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2723 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2724 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2725 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2726 if(.not.read2sigma.and.start_from_model) then
2727 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2728 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2729 start_from_model=.false.
2731 if(start_from_model .and. (me.eq.king .or. .not. out1file))
2732 & write(iout,*) 'START_FROM_MODELS is ON'
2733 if(start_from_model .and. rest) then
2734 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2735 write(iout,*) 'START_FROM_MODELS is OFF'
2736 write(iout,*) 'remove restart keyword from input'
2739 if (homol_nset.gt.1)then
2740 call card_concat(controlcard)
2741 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2742 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2743 write(iout,*) "iset homology_weight "
2745 write(iout,*) i,waga_homology(i)
2748 iset=mod(kolor,homol_nset)+1
2751 waga_homology(1)=1.0
2754 cd write (iout,*) "nnt",nnt," nct",nct
2761 c write(iout,*) 'nnt=',nnt,'nct=',nct
2764 do k=1,constr_homology
2777 do k=1,constr_homology
2779 read(inp,'(a)') pdbfile
2780 if(me.eq.king .or. .not. out1file)
2781 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2782 & pdbfile(:ilen(pdbfile))
2783 open(ipdbin,file=pdbfile,status='old',err=33)
2785 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2786 & pdbfile(:ilen(pdbfile))
2789 c print *,'Begin reading pdb data'
2791 c Files containing res sim or local scores (former containing sigmas)
2794 write(kic2,'(bz,i2.2)') k
2796 tpl_k_rescore="template"//kic2//".sco"
2799 if (read2sigma) then
2800 call readpdb_template(k)
2805 c Distance restraints
2808 C Copy the coordinates from reference coordinates (?)
2812 c write (iout,*) "c(",j,i,") =",c(j,i)
2816 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2818 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2819 open (ientin,file=tpl_k_rescore,status='old')
2820 if (nnt.gt.1) rescore(k,1)=0.0d0
2821 do irec=nnt,nct ! loop for reading res sim
2822 if (read2sigma) then
2823 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2824 & rescore3_tmp,idomain_tmp
2826 idomain(k,i_tmp)=idomain_tmp
2827 rescore(k,i_tmp)=rescore_tmp
2828 rescore2(k,i_tmp)=rescore2_tmp
2829 rescore3(k,i_tmp)=rescore3_tmp
2830 if (.not. out1file .or. me.eq.king)
2831 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
2832 & i_tmp,rescore2_tmp,rescore_tmp,
2833 & rescore3_tmp,idomain_tmp
2836 read (ientin,*,end=1401) rescore_tmp
2838 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2839 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2840 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2845 if (waga_dist.ne.0.0d0) then
2853 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2854 c write (iout,*) k,i,j,distal,dist2_cut
2856 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2857 & .and. distal.le.dist2_cut ) then
2863 c write (iout,*) "k",k
2864 c write (iout,*) "i",i," j",j," constr_homology",
2869 if (read2sigma) then
2872 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2874 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2875 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2876 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2878 if (odl(k,ii).le.dist_cut) then
2879 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2882 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2883 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2885 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2886 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2890 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2893 l_homo(k,ii)=.false.
2900 c Theta, dihedral and SC retraints
2902 if (waga_angle.gt.0.0d0) then
2903 c open (ientin,file=tpl_k_sigma_dih,status='old')
2904 c do irec=1,maxres-3 ! loop for reading sigma_dih
2905 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2906 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2907 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2908 c & sigma_dih(k,i+nnt-1)
2913 if (idomain(k,i).eq.0) then
2917 dih(k,i)=phiref(i) ! right?
2918 c read (ientin,*) sigma_dih(k,i) ! original variant
2919 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2920 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2921 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2922 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2923 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2925 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2926 & rescore(k,i-2)+rescore(k,i-3))/4.0
2927 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2928 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2929 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2930 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2931 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2932 c Instead of res sim other local measure of b/b str reliability possible
2933 if (sigma_dih(k,i).ne.0)
2934 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2935 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2940 if (waga_theta.gt.0.0d0) then
2941 c open (ientin,file=tpl_k_sigma_theta,status='old')
2942 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2943 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2944 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2945 c & sigma_theta(k,i+nnt-1)
2950 do i = nnt+2,nct ! right? without parallel.
2951 c do i = i=1,nres ! alternative for bounds acc to readpdb?
2952 c do i=ithet_start,ithet_end ! with FG parallel.
2953 if (idomain(k,i).eq.0) then
2954 sigma_theta(k,i)=0.0
2957 thetatpl(k,i)=thetaref(i)
2958 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2959 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2960 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2961 c & "rescore(",k,i-2,") =",rescore(k,i-2)
2962 c read (ientin,*) sigma_theta(k,i) ! 1st variant
2963 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
2964 & rescore(k,i-2))/3.0
2965 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2966 if (sigma_theta(k,i).ne.0)
2967 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2969 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2970 c rescore(k,i-2) ! right expression ?
2971 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2975 if (waga_d.gt.0.0d0) then
2976 c open (ientin,file=tpl_k_sigma_d,status='old')
2977 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2978 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2979 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2980 c & sigma_d(k,i+nnt-1)
2984 do i = nnt,nct ! right? without parallel.
2985 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
2986 c do i=loc_start,loc_end ! with FG parallel.
2987 if (itype(i).eq.10) cycle
2988 if (idomain(k,i).eq.0 ) then
2995 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2996 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2997 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2998 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
2999 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3000 if (sigma_d(k,i).ne.0)
3001 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3003 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3004 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3005 c read (ientin,*) sigma_d(k,i) ! 1st variant
3010 c remove distance restraints not used in any model from the list
3011 c shift data in all arrays
3013 if (waga_dist.ne.0.0d0) then
3019 if (ii_in_use(ii).eq.0.and.liiflag) then
3023 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3024 & .not.liiflag.and.ii.eq.lim_odl) then
3025 if (ii.eq.lim_odl) then
3026 iishift=ii-iistart+1
3031 do ki=iistart,lim_odl-iishift
3032 ires_homo(ki)=ires_homo(ki+iishift)
3033 jres_homo(ki)=jres_homo(ki+iishift)
3034 ii_in_use(ki)=ii_in_use(ki+iishift)
3035 do k=1,constr_homology
3036 odl(k,ki)=odl(k,ki+iishift)
3037 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3038 l_homo(k,ki)=l_homo(k,ki+iishift)
3042 lim_odl=lim_odl-iishift
3047 if (constr_homology.gt.0) call homology_partition
3048 if (constr_homology.gt.0) call init_int_table
3049 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3050 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3054 if (.not.lprn) return
3055 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3056 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3057 write (iout,*) "Distance restraints from templates"
3059 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3060 & ii,ires_homo(ii),jres_homo(ii),
3061 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3062 & ki=1,constr_homology)
3064 write (iout,*) "Dihedral angle restraints from templates"
3066 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3067 & (rad2deg*dih(ki,i),
3068 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3070 write (iout,*) "Virtual-bond angle restraints from templates"
3072 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3073 & (rad2deg*thetatpl(ki,i),
3074 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3076 write (iout,*) "SC restraints from templates"
3078 write(iout,'(i5,100(4f8.2,4x))') i,
3079 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3080 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3083 c -----------------------------------------------------------------
3086 c----------------------------------------------------------------------
3089 subroutine flush(iu)
3094 subroutine flush(iu)
3099 c------------------------------------------------------------------------------
3100 subroutine copy_to_tmp(source)
3101 include "DIMENSIONS"
3102 include "COMMON.IOUNITS"
3103 character*(*) source
3104 character* 256 tmpfile
3108 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3109 inquire(file=tmpfile,exist=ex)
3111 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3112 & " to temporary directory..."
3113 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3114 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3118 c------------------------------------------------------------------------------
3119 subroutine move_from_tmp(source)
3120 include "DIMENSIONS"
3121 include "COMMON.IOUNITS"
3122 character*(*) source
3125 write (*,*) "Moving ",source(:ilen(source)),
3126 & " from temporary directory to working directory"
3127 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3128 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3131 c------------------------------------------------------------------------------
3132 subroutine random_init(seed)
3134 C Initialize random number generator
3136 implicit real*8 (a-h,o-z)
3137 include 'DIMENSIONS'
3143 logical OKRandom, prng_restart
3145 integer iseed_array(4)
3147 include 'COMMON.IOUNITS'
3148 include 'COMMON.TIME1'
3149 include 'COMMON.THREAD'
3150 include 'COMMON.SBRIDGE'
3151 include 'COMMON.CONTROL'
3152 include 'COMMON.MCM'
3153 include 'COMMON.MAP'
3154 include 'COMMON.HEADER'
3155 include 'COMMON.CSA'
3156 include 'COMMON.CHAIN'
3157 include 'COMMON.MUCA'
3159 include 'COMMON.FFIELD'
3160 include 'COMMON.SETUP'
3161 iseed=-dint(dabs(seed))
3162 if (iseed.eq.0) then
3163 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3164 & 'Random seed undefined. The program will stop.'
3165 write (*,'(/80(1h*)/20x,a/80(1h*))')
3166 & 'Random seed undefined. The program will stop.'
3168 call mpi_finalize(mpi_comm_world,ierr)
3170 stop 'Bad random seed.'
3173 if (fg_rank.eq.0) then
3177 if(me.eq.king .or. .not. out1file)
3178 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3179 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3180 OKRandom = prng_restart(me,iseedi8)
3183 tmp=65536.0d0**(4-i)
3184 iseed_array(i) = dint(seed/tmp)
3185 seed=seed-iseed_array(i)*tmp
3188 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3189 & (iseed_array(i),i=1,4)
3190 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3191 & (iseed_array(i),i=1,4)
3192 OKRandom = prng_restart(me,iseed_array)
3195 c r1 = prng_next(me)
3196 r1=ran_number(0.0D0,1.0D0)
3197 if(me.eq.king .or. .not. out1file)
3198 & write (iout,*) 'ran_num',r1
3199 if (r1.lt.0.0d0) OKRandom=.false.
3201 if (.not.OKRandom) then
3202 write (iout,*) 'PRNG IS NOT WORKING!!!'
3203 print *,'PRNG IS NOT WORKING!!!'
3206 call mpi_abort(mpi_comm_world,error_msg,ierr)
3209 write (iout,*) 'too many processors for parallel prng'
3210 write (*,*) 'too many processors for parallel prng'
3218 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)